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ABSTRACT 


A  study  is  made  on  applying  the  turbulent  kinetic  energy  approach 
to  inhomogeneous  two-stream  turbulent  mixing  calculations  and  to  calcu¬ 
lations  of  a  two-dimensional  symmetric  wake.  Mixing  calculations  are 
made  and  compared  with  experimental  data  for  coaxial  hydrogen-air  and 
air-air  mixing  flows.  The  turbulent  kinetic  energy  equation  is  trans¬ 
formed  into  a  transport  equation  for  the  turbulent  shear  stress  by  assum¬ 
ing  that  the  turbulent  shear  stress  is  directly  proportional  to  the  turbulent 
kinetic  energy.  A  flux  model  is  assumed  for  the  lateral  diffusion  of  tur¬ 
bulent  kinetic  energy.  The  energy  dissipation  is  modeled  to  a  form  simi¬ 
lar  to  that  derived  for  isotropic  turbulence.  Mass  and  energy  transport 
are  incorporated  in  the  analysis  by  assuming  the  Prandtl  and  Lewis  num¬ 
bers  to  be  unity.  The  resulting  set  of  partial  differential  equations  is 
hyperbolic  and  the  method  of  characteristics  was  chosen  for  their  solu¬ 
tion.  The  theoretical  method  inherently  incorporates  history  of  the  tur¬ 
bulent  structure  in  the  calculations  which  is  physically  more  acceptable 
than  turbulent  structure  models  based  on  local  flow  properties.  The 
results  show  that  the  turbulent  kinetic  energy  approach  is  quite  applic¬ 
able  to  two-stream  mixing  problems,  although  the  method  requires  fur¬ 
ther  development  before  it  is  useful  for  routine  engineering  calculations. 
Calculations  for  a  wake  behind  a  flat  plate  are  made  and  found  to  compare 
well  with  experimental  data.  The  flux  model  for  the  diffusion  of  turbulent 
kinetic  energy  and  the  energy  dissipation  model  were  found  to  produce 
results  which  agree  well  with  experimental  data  for  both  the  mixing  jet 
and  the  wake. 
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CHAPTER  I 

I 

INTRODUCTION 

Turbulent  mixing  in  jets  and  wakes  has  been  studied 
for  several  decades.  Interest  in  inhomogeneous  two-stream 
turbulent  mixing  has  increased  significantly  in  the  last 
several  years  as  the  result  of  increased  interest  in  super¬ 
sonic  combustion,  gas  ejectors,  and  exhaust  jets.  At  best 
only  a  qualitative  understanding  of  the  physical  processes 
occurring  in  the  most  simple  turbulent  shear  flows  is  exis¬ 
tent.  Both  analytical  and  experimental  progress  are  still 
not  satisfactory. 

Until  recently  the  analytical  approaches  on  turbulent 
jet  mixing  could  be  classified  [1]^  as  follows: 

*  1.  Point  source  diffusion  of  momentum,  material,  or 

temperature  using  equations  and  solutions  well  known  from 
the  study  of  heat  flow. 

2.  Boundary- layer  form  of  the  Navier-Stokes  equation, 
into  which  are  inserted  various  transport  theories  such  as 

a.  Momentum  transport,  using  the  mixing  length 
concept 

b.  Vorticity  transport,  using  the  mixing  length 
concept 

c.  Constant  exchange  coefficient 

^Numbers  in  brackets  refer  to  similarly  numbered 
references  in  the  bibliography. 
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d.  Karman  similarity  theory 

e.  Statistical  theory 

3.  Integral  equations  of  motion. 

Also ,  Forstall  and  Shapiro  [1]  have  compiled  a  very  extensive 
list  of  references  covering  experimental  and  analytical  work 
on  jet  mixing  prior  to  1950.  Some  of  the  more  recent  work 
is  reported  in  References  2  through  21.  Most  analytical 
approaches  have  dealt  with  method  2  choosing  various  trans¬ 
port  theories.  Perhaps  the  most  widely  used  and  successful 

i 

transport  theories  are  those  proposed  by  Prandtl,  i.e.,  his 
mixing  length  theory  and  constant  exchange  coefficient 
hypothesis  discussed  in  Reference  22.  Both  the  mixing 
length  theory  and  constant  exchange  coefficient  hypothesis 
relate  the  structure  of  turbulence  through  the  shear  stress 
to  local  mean  flow  properties  which  in  effect  relate  the 
turbulence  structure  to  local  flow  properties.  These 
approaches,  however,  are  used  successfully  in  many  jet  mix¬ 
ing  problems.  Since  they  do  not  include  the  physics  of  past 
history  of  the  turbulent  structure,  they  do  not  account  for 

the  characteristics  of  the  experimental  apparatus  which 

/ 

induces  turbulent  history  into  the  flow.  Ferri  [23]  dis¬ 
cusses  this  problem  in  detail  particularly  as  related  to  the 
two-stream  mixing  problem.  The  characteristics  of  the 
experimental  apparatus  upstream  of  the  mixing  region 
(initial  boundary  layers)  have  important  effects  on  the 
mixing  phenomena  near  the  origin  of  mixing.  Oftentimes 
this  point  is  neglected  in  presenting  and  comparing 
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experimental  results.  Thus,  because  the  character  of  the 
initial  conditions  does  significantly  contribute  to  the 
mixing  phenomena  and  turbulence  structure  in  the  vicinity  of 
the  origin  of  mixing,  it  seems  reasonable  that  the  analytical 
solution  of  the  jet  mixing  problem  should  be  an  initial 
valued  problem  for  both  the  turbulent  and  mean  flow  fields. 
Such  an  analytical  approach  will  be  discussed  in  detail  in 
this  dissertation. 

The  basic  idea  was  derived  from  Bradshaw,  Ferriss,  and 
Atwell's  analytical  treatment  of  the  equilibrium  turbulent 
boundary- layer  problem  (24  and  25].  They  used  Townsend's 
suggestion  [26]  that  the  turbulent  shear  stress  was  propor¬ 
tional  to  the  turbulent  kinetic  energy  and  then  rewrote  the 
turbulent  kinetic  energy  equation  as  a  transport  equation 
for  the  turbulent  shear  stress.  They  also  assumed  that  the 
diffusion  of  turbulence  was  purely  convective  and  dominated 
by  the  large-scale  flow  eddies  as  opposed  to  small-scale 
h'igh-intensity  turbulent  motion.  Customarily  [27]  it  is 
assumed  that  the  small-scale  high— intensity  turbulence  dif¬ 
fuses  according  to  a  gradient  model;  i.e.. 


3c?2 

Diffusion  ~  -k-4— 
3y 


and  large-scale  eddies  diffuse  according  to 


Diffusion  ~ 

m  '  '  ***  N 

where  is  the  effective  rate  at  which  the  turbulence 
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diffuses  or  is  convected  in  the  lateral  direction.  Bradshaw 
et  al.  [28]  further  emphasizes  that  large  eddies  entirely 
dominate  the  free  mixing  layer.  If  this  is  correct,  a  con¬ 
vective  model  for  the  lateral  diffusion  appears  more  suited 
to  the  free- jet  mixing  problem  than  a  gradient  model. 

However,  Patankar  and  Spaulding  [29]  have  formulated 
the  problem  by  assuming  a  gradient  model  for  the  diffusion 
of  turbulent  kinetic  energy.  Along  with  this  assumption 
they  have  also  assumed  a  gradient  model  for  the  turbulent 
shear  stress.  These  two  assumptions  result  in  the  momentum 
and  turbulent  kinetic  energy  equations  both  being  parabolic. 
Lee  and  Harsha  [30]  have  studied  the  two-dimensional  wake 
and  the  two-stream  mixing  problems  utilizing  Bradshaw's 
assumption  [24]  that  the  turbulent  shear  stress  is  linearly 
proportional  to  the  turbulent  kinetic  energy  and  a  gradient 
model  for  the  diffusion  of  turbulent  kinetic  energy.  They 
used  the  Patankar  and  Spaulding  method  for  solving  the 
resulting  general  set  of  parabolic  equations.  The  method  is 
reasonably  successful.  In  the  Bradshaw  et  al.  approach  [24] 
the  shear  stress  is  not  modeled  in  the  momentum  equation, 
but  left  as  a  dependent  variable.  With  the  convective 
lateral  diffusion  formulation,  the  continuity,  momentum,  and 
turbulent  kinetic  energy  equations  form  a  first-order  quasi- 
linear  hyperbolic  set  of  partial  differential  equations.  In 
fact,  the  equations  are  still  hyperbolic  if  diffusion  is 
neglected.  Using  either  the  Lee  and  Harsha  approach  [30]  or 
Bradshaw's  approach  [24]  with  the  turbulent  kinetic  energy 
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equation  achieves  the  objective  of  formulating  the  jet 

mixing  problem  as  an  initial  valued  problem  with  respect  to 

- — — — *  ~ 

the  turbulent  structure .  /The  approach  in  this  dissertation 
is  to  use  the  Bradshaw  method  for  the  two-stream  mixing 
problem  and  incorporate  the  density  variation  by  use  of  the^ 


-Giroc.aoL^law^^jromparisons  with  air-air  and  hydrogen-air 
coaxial  mixing  experimental  data  will  be  made.  A  brief 


study  is  also  presented  on  the  use  of  the  Bradshaw  method  to 
calculate  wake  flow  behind  a  thin  flat  plate. 
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CHAPTER  II 

DERIVATION  OF  THE  COMPRESSIBLE  TURBULENT 
KINETIC  ENERGY  EQUATION 

I.  TIME  AVERAGE  RULES 

A  derivation  of  the  compressible  turbulent  kinetic 
energy  equation  will  be  made  along  with  legitimate  approxi¬ 
mations  to  simplify  the  equation  to  useful  form.  The  equa¬ 
tion  can  be  developed  entirely  from  the  instantaneous  form 
of  the  continuity  and  the  momentum  equations  along  with  an 
averaging  procedure.  The  usual  Eulerian  time-averaged 
definition  is  used  [27].  The  time-averaged  value  of  any 
dependent  variable  is  defined  by 

t 

A  =  lim  f  A  dt 

t-  2t 

Thus  designation  of  the  overscored  quantity  as  a  time- 
averaged  quantity,  the  instantaneous  quantity  can  be  written 
as 

A  =  A  +  A' 

where  A'  is  the  instantaneous  fluctuation  away  from  the 
mean.  The  usual  rules  of  time  averaging  are  employed  for 
the  products  of  fluctuating  quantities.  These  are 

A  =  A  +  A'  =  A  +  A'  =  A  +  A' 
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thus 


A'  =  0 

With 

B  =  B  +  B' 

A  B'  =  A  B 1  =  A  B 1  =0  since  B'  =  0 

and 

A  B  =  (A  +  A'  }  (B  +  B '  > 

=  A  B  +  A  B'  +  B  A*  +  A1 B' 

=  A  B  +  A 1 B 1 

At  this  point  the  compressible  form  of  the  instan¬ 
taneous  equations  of  motion  will  be  stated  followed  by  their 
conversion  to  their  time-averaged  form.  By  manipulation  of 
the  instantaneous  and  time-averaged  equations,  the  compress¬ 
ible  form  of  the  turbulent  kinetic  energy  equation  will  be 
derived.  This  treatment  of  the  equation  will  be  accomplished 
in  Cartesian  coordinates  for  simplicity.  The  generalized 
two-dimensional  form  of  the  final  equations  is  used  in 
succeeding  chapters  for  both  the  plane  two-dimensional  wake 
and  the  axisymmetric  two-stream  jet  mixing  problems. 

II.  CONTINUITY 

The  instantaneous  continuity  equation  is 


7 


AEDC -TR-70-134 


It  +  377  (pui}  =  0 


(1) 


Using  the  time-averaged  definition 


P  =  P  +  P' 


U.  =  U.  +  U! 
ill 


where  (')  is  used  to  designate  the  fluctuating  quantity. 
Substituting  into  Equation  1  and  time  averaging  gives 


sir  <p  Di  +  <,'Di)  "  0 


(2) 


which  is  the  steady-state  time-averaged  form  of  the  con¬ 
tinuity  equation. 

Ill .  MOMENTUM 

The  instantaneous  momentum  equation  is 


3U,  3U . 


3P 

8xi 


3x  i 
D 


(3) 


where  is  the  viscous  stress  tensor.  Substituting  in  the 
instantaneous  fluctuating  properties  we  have 


3  <U.  +  U! ) 


_  3(U.  +  U[) 

(P  +  p«) - .  1  +  (p  +  p')  (Uj  +  up — X3X| 


3P  3P 1 


ax. 


3x^  •.  3Xj 


3cjj_  _  ffii 


3x . 
3 


(4) 


8 


AEDC-TR-70-134 


Taking  the  time  average  gives 


_  3u .  _  3u !  ~  3ur 

(p  uj  +  p'uj>  inr  +  c  uj  h1  +  »'  (uj  -  “jhnr 

■1  3  3 


3P  3ali 


3x^  3Xj 


(5) 


By  multiplying  the  instantaneous  continuity  equation  by  the 
instantaneous  velocity,  U^#  time  averaging  and  using 
Equation  2  gives 


ui  sir  (p  Uj  +  p'Uj  +  =  0 


(6) 


By  adding  this  term  to  Equation  5  and  combining  terms 


( 


3U .  dp  U'.U'.  Sp'u^.u:  3p’u!  u. 

p  uj  +  5^  +  +  +  -  ax,  3 


ap_  -  !fii 

3x.  3x . 

i  1 


(7) 


which  is  the  complete  steady-state  time-averaged  form  of  the 
momentum  equation. 

IV.  TURBULENT  KINETIC  ENERGY 


First  the  instantaneous  momentum  equation  is  multi¬ 
plied  by  U^  and  instantaneous  continuity  equation  by  1/2  U?. 
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Continuity 


Uf 


U? 


"i  It +  “i  db  (pUj} 


=  0 


Momentum 


3U?  pU.  3U? 

i  “  -i  * 


jp  _ i 

2  3t 


2  3x. 

3 


-  u.  |P_  .  D.  !!ii 

l  3x^  l  3Xj 


Adding  these  two  equations  together  gives 


3 

5t 


'pu?' 


2  3x 


.  u  .  u  !!ii 

ui  8x±  ui  8x. 


Time  averaging  the  steady-state  form  of  this  equation  leaves 


1  3  “Vi  7- n_ 


3Xj 


5^-7 

_  n  — _ u  ±J- 

ui  3x±  ui  3Xj 


(8) 


Expanding  the  terms  of  this  equation 


I  sir  tp<uiuj  +  VT  +  2  ui“Fj  +  iTuj)1 


1  3 


+  j  sir  ‘2  uiuj  +  uj  +  UI 


+  2  u.  p'Uiu’.  +  p ■  u :  !u ! ] 
1  1  j  1  j 


do.  . 


do]T 


U  .  v  Ml  -  u  _Ai  -  u-  -ii 

ui  3x±  ui  3xi  ui  3x_.  ui  3x.. 


(8a) 


Multiply  Equation  7  by  U.  and  subtract  from  Equation  8a 
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U?  dp  U.  U?  3p'u:  3U . 

_L_1  +  -1—1  +  p  uTDT  ^ 

2  3x.  2  3x.  v  l  j  3x . 

3  3  3 


3U .  3U .  ,  3p  U .U! ‘ 

+  uj  asr  +  2  ~  iif- 


,  3p  U'2U!  ,  a _ ,  a  - 

+  4  - r-i - 1  +  4  ——  p  1 U  !  2U  .  +  4  7~-  p  'U  !  2U 

2  3x.  2  3x.  F  l  T  2  3x.  w  l  3 

3  3  3 


3  a 


=  _  a»  «1.  D!  pi 

1  3Xj  1  3Xj 


(9) 


The  first  two  terms  are  identically  zero  as  given  by  the 
time-averaged  continuity  equation.  Changing  the  suffix  on 
this  pressure  term  from  i  to  j  ,  and  after  some  rearranging, 
the  compressible  form  of  the  turbulent  kinetic  energy 
equation  is 


a  f-i  , _  3U. 

uj  (W2  +  uj  55E7 


II 


.  . _  ^30.  3U 

|oUi2  +  ^ 


III 


3 

IV 


3U. 


+  (p  tiTuT  +  pW)  ^ 


3x 


-  +  +  ip'ui'uj]  +  Di  -sj- = 


=  0  (10) 


VI 


VII 
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The  bars  over  the  mean  flow  quantities  have  been  dropped  for 
convenience.  The  seven  labeled  terms- of  the  turbulent 
kinetic  energy  are: 

I .  Convection  by  mean  flow 

II.  Turbulent  mass  flux  times  mean  flow  acceleration 

III.  Normal  stress  times  mean  dilatation 

IV.  Fluctuating  pressure  dilatation 

V.  Turbulence  production 

VI.  Turbulent  energy  diffusion 

VII.  Viscous  dissipation  of  turbulent  energy 
Equation  10  is  an  exact  expression  for  the  turbulent  kinetic 
energy  equation  including  compressibility.  An  order  of 


magnitude  analysis  will  be  performed  and  simplifying 
approximations  made  to  reduce  this  equation,  and  the  con 


^tinuity  andmomentum  equations,  to  a  more  usajble—  form.._ _ 

First,  an  approximation  to  the  fluctuating  density,  p',  is  J 
developed.  ^The  instantaneous  total  enthalpy  and  equation 


of  state  for  aV^thermally  and  calorically  perfect  gas__areJ 
given  by 


The  instantaneous  properties  are  defined 
H  =  H  +  H' 
h  -  h  +  h’ 
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U  =  U  +  U' 

P  =  P  +  P' 

p  =  p  +  p* 

With  these  definitions ,  Equation  11  is  written  as 
H  +  H1  =  h  +  h'  (U?  +  2U.U!  +  U!2) 

»  X  JL  X  JL 

Time  averaging  and  taking  the  difference  gives 

U! 2  -  U! 2 

H'  =  h'  +  U.U!  +  -i - = — — 

li  2 

The  last  term  is  assumed  to  be  an  order  of  magnitude  less 
than  U.U!.  Also  U' ,  V',  and  W'  are  assumed  to  be  of  equal 
order  and 


h'  +  U  U'  -  0  (13) 

From  the  equation  of  state 


P  +  P'  =  [p  H"  +  p'E  +  ph'  +  p'h']  (14) 

The  time  average  of  Equation  14  is 

P  =  [p  h  +  VlF]  (15) 

Subtracting  Equation  15  from  Equation  14  leaves 
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P'  =  ( P ' h  +  ph'  +  p'h'  -  pF7*) 


Bradshaw  [25]  suggests  that  the  first  and  second 
terms  on  the  right-hand  side  are  the  dominant  terms 
when  the  Mach  number  fluctuation  is  much  less  than  unity. 
This  gives 


£l  a  _  hi 

p"  h 


(16) 


Eliminating  h1 'betwelen  Equations  13  and  16, 


£l  ~  0  U1 

P  h 


or 


,  (y  -  i)  M2  ^ 
p  U 


where_M.is.the  local  time-averaged  Mach  number.  This 


\ 


(17)  ! 

v  /' 


J 


(approximation  will  be  used  to  simplify  the  contini  /ty,  . 
momentum,  and  turbulent  kinetic  energy  equations. /  It  is 
also  assumed  that  (y  -  1)M2  is  no  greater  than  order  unity. 
The  reader  is  referred  to  Appendix  A  for  the  order  of  magni¬ 
tude  analysis  of  the  continuity,  momentum,  and  turbulent 
kinetic  energy  equations.  The  resulting  equations  are: 

Continuity 


|£U  +  |_  tpv  + 


=  0 


(18) 
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3U 


pU  +  (pV  +  P’V1) 


TTrT '  3U 


37 


3P  3 


3x  3y 


(pU'V'  +  p'U'V' ) 


(19) 


Turbulent  Kinetic  Energy 


2  ,^T  j.  :rr?m  3^2  ,  -^rrrrrK  9U 


P£  Ml  +  (PV  +  P'V  )  0S_  x  /onTTTT  x  7H nT,TT 

/\  rv  1  n 


2  3x 


3y 


+  (plJrVr  +  p’u'V1  ) 


3y 


+  l^r  (FF  +  \  pq^'  +  |  p  V2V<)  +  u|  ^±1 


TaH 


3y 


(20) 


In  a  manner  similar  to  Bradshaw  [24]  the  following 
definitions  are  made: 


al  = 


P?2 


(21) 


Vd  E 


(P'V1  +  i-  pq^^V'  +  \  p  ' q '  *V') 


pq' 2/2 


(22) 


r  "i 


L  = 


axp 


3/2 


(23) 


3a!  . 

U!  *1 
l  3x . 
3 


Using  these  definitions,  neglecting  the  viscous  stresses, 
letting  the  variable  V  +  p 1 V'/p  he  defined  as  V,  and  defin 
ing  the  turbulent  shear  stress  as 
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T  =  -  (pU'V'  +  p'U'V' ) 


the  continuity,  momentum,  and  turbulent  kinetic  energy 
equations  become 


i ££  +  lev 

3x  3y  u 


(24) 


..  3U  ,  ..  3U  ... 

pU  33?  +  PV  37  " 


3P  3t 
3x  3y 


(25) 


3 

PU  _ 

2  3x 


(-E-1  ,U.l 
Uxpj  pV  laxpj 


3y 


-  T 


3U 


3y 


Vd  T 

2~  FT 


■J 


f— 1 

laxPj 


3/2 


+  pD  =  0 


(26) 


The  last  term,  pD,  in  Equation  26  is  the  same  as  the  last 
two  terms  in  Equation  20.  Using  the  approximations  given  by 
Equation  17,  the  function  D  can  be  estimated  as  follows: 


P^U 


T 


pU'  2 


u 

h 


fFV1*  = 


tU 

h 


It  has  been  noted  by  References  21  and  26  that  U'2  cJ2/2. 

Then 


p'U' 


T  U 
2a^  h 


With  these  relations  and  by  use  of  the  momentum  equation,  pD 
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is  approximated  as 


T  Uf 

-  + 

111 

9 

rxU 

T 

2a^p  h[ 

9x 

3yJ 

"  3y 

hj 

[2a1pJ 

This  expression  is  rearranged  to  a  more  usable  form: 


pD 


_ r_  U  9P 

2a^p  h  9x 


u  a_ 

T  1 

h  9y 

i 

[  2a1pJ 

+ 


T  9 _ fu  j 

2axp  9y (hj 


(27) 


(28) 


This  completes  the  development  of  the  basic  equations 
of  motion  in  two-dimensional  flow.  The  subsequent  analysis 
in  Chapter  III  uses  the  generalized  two-dimensional 
equations  (i.e.,  including  axisymmetric  flow). 
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CHAPTER  III 

CHARACTERISTIC  AND  COMPATIBILITY  EQUATIONS 

I.  BASIC  ASSUMPTIONS  AND  DEVELOPMENT 


If  real- characteristics  exist  for  Equations  24,  25, 
and  26,  then  this  system  of  equations  is  hyperbolic  [32]. 
However,  prior  to  investigating  the  characteristics,  addi¬ 
tional  assumptions  will  be  made  to  facilitate  the  develop¬ 
ment  of  the  equations: 

a.  The  flow  is  plane  two-dimensional  or  axisymmetric . 

b.  The  mixing  process  is  a  parallel  two-stream 
turbulent  mixing  process. 

c.  The  turbulent  Prandtl  and  Schmidt  numbers  are 
unity,  implying  that  the  velocity,  energy,  and  mass  species 
profiles  are  similar.  Very  near  the  initial  plane  of  mixing 
where  the  mixing  zone  is  undeveloped  and  dominated  by  the 
initial  boundary  layers,  profile  similarity  does  not  exist. 
Thus,  this  assumption  more  properly  applies  to  the  fully 
developed  mixing  zone. 

Equations  24,  25,  and  26  are  restated  and  used  in 
generalized  two-dimensional  coordinates  as  follows: 


Mass  Continuity 


3pUrk  .  3pVrk 
3x  “Sr 


(29) 
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Momentum 


„  9U  .  „  3U 
pU  3?  +  pV  9?  ' 


dP  J.  3rkT 

dx  k  3r 

r 


(30) 


Turbulent  Kinetic  Energy 


r 

• 

f.T  3  ...  3  ,  _1  3_ 

[  3x  ■  3rJ  pa1  rk  3r 

k 

vd 

r 

2  ci  i 
1) 

T 

4 

-  T 


9U  p 
3r  L 


pa. 


3/2 


+  pD  =  0 


(31) 


where  k  is  zero  and  one  for  plane  two-dimensional  flow  and 
axisymmetric  flow,  respectively.  The  term  D  is  still  given 
by  Equation  28. 

If  U,  V,  and  r  are  chosen  as  the  dependent  variables 
in  Equations  29,  30,  and  31,  an  expression  for  the  density 
must  be  developed  in  terms  of  one  or  more  of  these  three 
dependent  variables.  The  density  can  be  expressed  in  terms 
of  U  by  use  of  the  Crocco  law  since  we  have  assumed  the 
Prandtl  and  Schmidt  numbers  to  be  unity.  From  the  equation 
of  state 


P 


(32) 


where  M  is  the  average  molecular  weight  and  Rq  is  the 
universal  gas  constant.  The  local  gas  mixture  is  assumed  to 
be  calorically  perfect,  then 


P 


=  llSe. 

*oh 


(33) 
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where  Cp  is  the  average  molar  specific  heat  and  h  is  the 
static  enthalpy.  Thus,  in  a  two-stream  compressible  mixing 
process  the  density  is  a  function  of  pressure,  the  molar 
specific  heat,  and  enthalpy.  The  molar  specific  heat  and 
enthalpy  can  be  expressed  in  terms  of  velocity  as  will  be 
seen.  The  axial  pressure  variation  must  be  given.  The  log 
differential  of  Equation  33  is 

1  3p_  1_  _  1  3h_  1  3P_  n4 

p  3x±  Cp  3x^  ”  h  3x.^  P  3x^  ' 

where  x^  represents  x  and  x2  represents  r.  The  average 
molar  specific  heat  is 


C  =  E  Y  C 
P  a  P 


(35) 


a 


where  Yq  is  the  mole  fraction  of  specie  a  given  by 


v  na 

Y«-  K 


(36) 


and  na  is  the  moles  of  specie  a.  For  two-stream  mixing 
where  only  two  species  need  be  considered 


n. 


n, 


*1  *  ? 

c  =  - ± -  c  +  — =■ —  c 

p  ni  +  n2  Pi  ni  +  n2  ^2 


Lettinc  C  be  the  mass  fraction 
a 


ci  _  niMi 
^2  ^2^2 


(37) 


(38) 
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and  substituting  into  Equation  37  gives 


C  = 
P 


Cn  M_C  +  C-M..C 

1  2  p1  2  1  p2 

C2M1  +  ClM2 


(39) 


Since  C1  +  C2  =  1 


C  = 
P 


Cl(M2CPl  '  "1CP2>  +  M1CP: 
ci(m2  -  Ml)  +  ' 


(40) 


Differentiating  Equation  40 


3C 

__E 

9xi 


(M_C  -  M.C  )  ~  C  (M_  -  M.) 
2  Pi  1  P2  P  2  1 


C1(M2  "  Ml}  +  M1 


3C1 

JxT 

i 


and  dividing  by  C 


n  3C 

1 _ E 

C  3x . 
P  '  i 


9C 

MlM2(cPl  '  V  557 


2  — 1 

{  C.  (M,C  -  M.C  )  +  M,C  K  C,  (M, 


'1  '“2 


Pi  1  P2  1  P2  1  2 


M.^  +  M^} 
(41) 


Note  that  mixing  of  perfect  monatomic  or  diatomic  gases 
gives 


3C 


3x 


^  =  0 


since 


C  =  C 

Pi  P; 
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By  the  assumption  of  unity  turbulent  Schmidt  number  the 
Crocco  law  gives 


U 


U. 


(42) 


where  the  subscripts  o  and  j  refer  to  outer  and  inner  (jet) 
boundary  conditions,  respectively.  Differentiating  Equation 
42 


3^ 

fci  ■  -  C. 

It  lo 

3U 

3x7  - 

u.  -  u 

3x . 

T  O  . 

J  / 

i 

(43) 


For  convenience  define 


M.M-(C  -  C  ) 
i  2  Pl  p2 


cii  "  Cl° 


F,  = 


U  .  -  U 
3  o 


lCl(M2CPl 


-  )  +  M.C  ]  [C,  (M,  - 


2  P 


P2 


'1  2 


+  Mx) 


(44) 


Then  from  Equations  41, 


43,  and  44 


(45) 


From  the  assumption  of  unity  turbulent  Prandtl 
number,  the  Crocco  law  gives 


H 

H 


j 


u 


u. 


where  H  is  the  total  enthalpy  defined  by 


(46) 
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H  =  h  + 


Uz  +  V' 


Assuming 


V 


<<  1 


H  =  h  +  U2/2 


(47) 


Differentiating  Equations  46  and  47 


3H  Hi  Ho  3U 

3x.  U.  -  U  3x. 

1  j  o  1 


(48) 


and 


3H  =  3h_  3U 

3x^  3x^  3x^ 


(49) 


By  eliminating  3H/3xi  between  Equations  48  and  49 


3h 

3x . 
i 


H.  -  H 

_] _ 2. 

,u .  -  u 

l  3  o 


-  U 


3U 

3x. 

l 


(50) 


Again  for  convenience  define 


F  =  — 
*  2  h 


H.  -  H. 


,U  .  -  U 

l  7  © 


-  U 


(51) 


Then 


1  3h_  3U__ 

h  3x .  2  3x. 

l  i 


(52) 


Combining  Equations  34,  45,  and  52  gives 


23 


AEDC-TR-70-134 


1  3p_ 
P  3XJ, 


-  F2) 


3U 

9x . 
1 


1  3P_ 

F  3x . 

l 


(53) 


This  is  the  expression  which  will  be  used  to  express  density 
gradients  in  terras  of  velocity  gradients.  Again,  recall 
that  equal  to  zero  corresponds  to  mixing  of  thermally  and 
calorically  perfect  monatomic  or  diatomic  gases.  Incom¬ 
pressible  homogeneous  mixing  corresponds  to  F2  equal  to 
zero . 

With  Equation  53  and  the  equation  of  state  the 
density  gradient  can  be  removed  from  explicity  appearing  in 
the  continuity,  momentum,  and  turbulent  kinetic  energy 
equations . 

First  the  continuity  equation  is 


U  +u(Fi  .  F  )]|£  +  V(Fl  -  F2)f 


+  ^  fe(rkv)  +  I  a!  “  0 


(54) 


The  momentum  equation  is 


V  -  i(Fl  -  F2) 


3U 

3r 


„  _  <Y  ~  1)  h  dP  .  1  3r  "(t/p) 
y  P  dx  k  3r 


(55) 


and  the  turbulent  kinetic  energy  equation  becomes 


ifu  L.  +  V  U 
2[u  3x  v  3rJ 


pa. 


k  3r 


V. 


2a, 


24 


AEDC-TR-70-1 34 


with 


_V 

2a, 


T 

PJ 


(F  .  F  )  22  -  I  22. 

lll  c2'  3r  p  3r 


pa. 


3/2 


+  D  =  0 


D  - 


23^ 


Y  -  1 


U  dP 
P  dx 


1  T  U 

2  p  h 


pa. 


3r 


1  "  F2U 

1  T 

T 

2  P 

axp 

•  . 

h 

SU 

3r 


(56) 


(57) 


Note  that  the  parameters  T  and  p  always  appear 
together  in  Equations  54,  55,  and  56  as  r/p.  For  convenience 
x/p  will  be  written  in  the  future  as  t. 

At  this  point  we  assume  the  functions  a^,  V^,  and  L 
are  functions  only  of  the  independent  variables  x  and  y.  It 
should  be  remarked  that  a^,  V^,  and  L  also  could  be  functions 
of  U,  V,  and  t  and  that  the  resulting  equations  would  be 
hyperbolic  if  the  characteristics  were  real.  However,  they 
cannot  in  general  be  functions  of  the  derivatives  of  the 
dependent  variables,  reference  Courant  and  Friedrichs  [32). 

The  method  used  for  the  development  of  the  character¬ 
istic  directions  and  compatibility  equatibns  is  described  in 
Appendix  B.  The  method  was  shown  to  the  author  by  Hr.  Fran 
Loper  of  ARO,  Inc.  Other  more  involved  methods  exist  and 
are  described  in  References  32  and  33.  The  three 
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characteristic  directions  are 


dr 

dx 


00 


(58) 


and 


dr 

dx 


L,R 


V  + 


Vd 

T 


—  —  (F  - 
2  U1 


F2) 


tU 

2h 


U 


V, 


+  I.  (p  _  p  )  -  “I 

2  2  2‘  2hJ 


+  T 


2al  '  Vd<Fl  -  F2>  +  (1  '  F2U) 


1/2 


T  U 


(59) 


One  characteristic  direction  is  normal  to  the  X-axis,  the 

left  characteristic  direction  is  at  an  angle  above  the 

■ 

streamline,  and  the  right  characteristic  direction  is  at  an 
angle  below  the  streamline.  It  is  noted  that  the  inclined 
characteristic  directions  are  real  provided 


tU 

"  2Kj 


2 


+  t  2ax  -  Vd  (F 


"  F2) 


l  <1  - 


f2u) 


>  0 


(60) 


For  incompressible  homogeneous  mixing  this  expression 
reduces  to 


Vd 

— j  +  2a^r  >  0 
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Both  terms  are  always  positive,  thus  the  system  of  equations 
is  hyperbolic. 

The  compatibility  equation  along  the  vertical 
characteristic  is 


[1  ♦  UCfj  -  F2)l 

-  [V  -  TIP,  -  Fa)  -  UT(P1  -  Fa>*lH+ 


+ 


u 


2  -  1 


-  1 


h[l  +  U  (F. 


"  F2)J 


'  1 
P 


dP 

dx 


=  0 


(61) 


This  equation  can  also  be  derived  directly  from  the  con' 
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G  =  2SJ  ^  3F(vd  rK) 


+ 


3/2 


T 

2a, 


-  llu  dP 

]P  3x 


UT  3ai 


2a 


2  3x 


V  +  Vd 


tU' 

“  “IT 


(63) 


This  completes  the  development  of  the  characteristics  and 
their  compatibility  equations. 


II .  BOUNDARY  CONDITIONS 


Symmetric  Flow 

For  the  symmetric  two-dimensional  flow,  first-  and 
second-regime  mixing  exist  as  illustrated  in  Figure  1.  The 
first-regime  boundary  conditions  ares 

Outer  Boundary 

T=  0 

u  =  uo 

V  =  vo 

vd  -  vdQ  164 1 

Inner  Boundary 
T  =  0 

u  =  u. 

1 
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V  =  0 


(65) 


At  both  the  inner  and  outer  boundaries  the  slopes  of  the 
right  and  left  characteristics  reduce  to 


dr 

b 

V. 

v  +  T-  1 

vd- 

2~ 

dx 

T.  t> 

u 

(66) 


as  can  be  seen  from  Equation  59.  The  superscript  b  is  used 
to  denote  the  boundary.  More  specifically,  at  the  outer 
boundary  with  positive. 


dr 

3x 


bo 


V  +  v, 

o  a( 

U 


(67) 


and 


dr 

dx 


bo 

R 


(68) 


\ 

The  left  characteristic  coincides  with  the  outer  boundary  of 
the  mixing  layer.  The  right  characteristic  coincides  with 
the . streamline .  At  the  inner  boundary,  with  negative. 


dr 

as 


bi 

L 


V. 

r 


0 


(69) 


and 
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dr 

as 


bi 

R 


(70) 


The  left  characteristic  coincides  with  the  streamline  at  the 
inner  boundary  and  the  right  characteristic  coincides  with 
the  edge  of  the  boundary.  If  the  inner  potential  core  flow 
is  one  dimensional,  the  left  characteristic  on  the  inner 
boundary  has  zero  slope. 

At  both  the  inner  and  outer  boundaries  each  term  of 
the  compatibility  equations  vanishes  (Equation  62) .  There¬ 
fore,  these  equations  provide  no  information  at  the  math¬ 
ematical  boundaries.  Normally,  this  would  present  some 
difficulty  in  solving  the  equations.  This  difficulty  is 
avoided,  however,  because  the  outer  boundary  is  arbitrarily 
defined  to  occur  at  the  position  where  the  U  component  of 
velocity  is  1  per  cent  of  the  maximum  velocity  difference  of 

the  two  streams.  This  criterion  is  for  U./U  greater  than 

j  o 

unity.  Likewise,  the  position  of  the  inner  boundary  is 
defined  to  occur  at  99  per  cent  of  the  velocity  difference. 

At  both  of  these  defined  boundaries  a  finite  value  of  shear 
stress  occurs.  Consequently,  the  compatibility  equations 
are  not  singular  at  these  defined  boundaries. 

In  the  second  regime  the  outer  boundary  conditions 
are  the  same  as  in  the  first  regime,  and  the  outer  boundary 
is  similarly  defined.  The  inner  boundary  conditions,  or  line 
of  symmetry  boundary  conditions,  are 

V  =  0 
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T  =  0 


U  = 


U 


3U 

3r 


0 


(71) 


The  shear  stress,  vertical  component  of  velocity,  and  diffu¬ 
sion  velocity  are  zero  by  symmetry.  These  conditions  give^/ 
the  slopes  of  both  the  right  and  left  characteristics  to  be 
zero  on  the  line  of  symmetry.  The  compatibility  equations 
are  singular  and  vanish  on  the  line  of  symmetry. 

The  only  unknown  on  the  line  of  symmetry  or  center- 

line  is  U  .  Thus,  one  equation  is  needed  to  calculate  U 
c  c 

and  is  obtained  from  either  the  continuity  or  momentum 
equation.  Using  the  momentum  equation  (Equation  55)  and 
taking  the  limit  as  r  approaches  zero  gives 


3U 

°c  TT--  iri T1 1  r  S  +  11  +  k>  H  (72> 

By  prescribing  the  pressure  gradient  and  calculating  the 
shear  stress  gradient  from  the  value  of  t  and  r  calculated 
from  the  first  characteristic  off  the  line  of  symmetry,  U 

C 

can  be  calculated. 


Unsymmetric  Flow 

If  the  flow  is  not  symmetric,  as  is  the  case  of  mix¬ 
ing  two  semi-infinite  streams,  the  vertical  component  of 
velocity  at  either  of  the  two  boundaries  is  not  known  a 


32 


AEDC-TR-70-134 


priori.  In  this  case,  the  vertical  components  of  velocity 
at  the  two  bouhdaries  are  related  to  one  another.  For 
incompressible  flow,  Mills  [34]  gives 


p.U.V.  +  p  U  V  =  0 
i  x  Ko  o  o 


(73) 


The  U  component  of  velocity  is  either  known  or  calculated 
from  the  potential  flow. 


4 
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CHAPTER  IV 

i 

PARAMETER  MODELING  FOR  FREE  JET-MIXING 

The  purpose  of  this  chapter  is  to  establish  physical¬ 
ly  perceptive  models  for  the  parameters  a^,  V^#  and  L  for 
axisymmetric-free  jet  mixing.  Already  these  parameters  have 
been  assumed  to  be  at  most  functions  of  the  independent 
variables  in  the  development  of  the  characteristic  and  com¬ 
patibility  equations.  Experimental  data  contained  in  the 
literature  will  be  used  as  a  guide  where  possible.  In 
Chapter  VI  the  models  established  here  are  used  in  the  cal¬ 
culation  of  several  cases  of  axisymmetric  jet  mixing  and  the 
calculations  are  compared  with  experimental  data. 

I.  TURBULENT  KINETIC  ENERGY-SHEAR 
STRESS  PARAMETER 

For  the  equilibrium  boundary  layer  Bradshaw  et  al. 

[24]  observed  that  a^  could  be  assumed  a  constant  equal  to 
0.15.  This  value  was  used  to  obtain  reasonable  results  in 
the  calculation  of  skin  friction  coefficient  and  velocity 
and  shear  stress  profiles  for  incompressible  and  compressible 
adiabatic  flows  with  and  without  pressure  gradient.  The 
choice  of  a^  equal  to  0.15  was  based  on  the  turbulent 
property  measurements  of  Klebanoff  [35] .  For  a  two- 
dimensional  plane  wake  and  an  axisymmetric  jet  Lee  and 
Harsha  [30]  observed  an  average  value  of  a^  equal  to  0.1 
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from  several  sources  [36,  37,  38,  and  39].  Recently  they 
have  revised  their  estimate  to  a^  equal  to  0.15  [40].  From 
the  turbulent  property  data  of  Sami  et  al.  [21]  for  a  homo¬ 
geneous  axisymmetric  jet  exhausting  into  still  air.  Figure  2 
was  prepared  to  show  the  variation  of  a^  with  the  non¬ 
dimens  ionali  zed  independent  variable  (r  -  r^)/b.  It  is 
shown  that  ja^j  can  be  approximated  reasonably  well  as  a 
function  of  (r  -  r^)/b  based  on  these  data.  The  data  shown 
are  taken  in  both  the  first  and  second  regimes  (3  <  x/D  <10). 
Wygnanski  and  Fiedler  [41]  have  recently  reported  other 
turbulent  intensity  measurements  in  an  incompressible  and 
axisymmetric  jet.  Their  measurements  give  a^  m  approxi¬ 
mately  equal  to  0.12  for  50<x/D<90.  Sami's  data  give  |ai|m 
approximately  equal  to  0.17  for  3  < x/D  < 10.  On  the  center- 
line  in  the  second  regime  a1  equals  zero  since  the  shear 
stress  is  zero,  but  the  turbulent  kinetic  energy  is  a  finite 
value.  The  algebraic  sign  of  a1  is  always  the  same  sign  as 
the  shear  stress  since  the  turbulent  kinetic  energy  is 
always  positive.  The  experimental  data  used  in  this  research 
for  comparison  with  the  theory  are  limited  to  x/D  <18;  thus, 
a  value  of  a^  based  on  Sami's  data  is  used  in  the  calcula¬ 
tions.  An  analytic  function  of  the  form 


Z  < 

o 


z  >  z 

O 


(74) 
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where 


r  -  r. 


Z  = 


is  used  in  the  calculations.  With 


ll  m 


equal  to  0.17  and 


Zq  equal  to  0.4  this  function  is  plotted  in  Figure  2  and  is 
seen  to  fit  Sami's  data  very  well  near  the  line  of  symmetry, 
In  Chapter  VI  results  of  the  calculations  with  other  values 
of  la, I  and  Z  are  discussed.  For  Z  >  Z  ,  a,  equal  to  a 
constant  was  found  to  be  satisfactory  in  all  calculations. 


II.  DIFFUSION  VELOCITY  PARAMETER 


The  parameter  V.  is  in  effect  the  rate  of  diffusion 
of  turbulent  kinetic  energy  in  the  lateral  direction.  By 
assuming  Vd  to  be  a  function  of  the  independent  variables, 
x  and  r,  the  continuity,  momentum,  and  turbulent  kinetic 
energy  equations  were  observed  in  Chapter  III  to  form  a 
hyperbolic  set  of  equations  for  the  three  dependent  variables 
U,  V,  and  t.  The  gradient  diffusion  models  create  a  para¬ 
bolic  set  of  partial  differential  equations.  The  justifies-: 
tion  for  the  nongradient  or  diffusion  velocity  concept  is 
based  on  the  idea  that  the  turbulent  diffusion  is  primarily 
dominated  by  the  effect  of  large  eddies  in  free  shear 
flow  [Z-®]  . 

a 

The  functional  form  of  the  parameter  V^  was  empirical¬ 
ly  derived  by  Bradshaw  et  al.  [24]  for  the  equilibrium 
boundary-layer  case  and  adjustments  were  made  by  trial  and 
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38 


error  until  plausible  solutions  were  obtained.  In  this 
study  a  similar  approach  is  taken  for  the  two-stream  mixing 

i 

problem.  The  value  of  can  be  empirically  observed  at  the 
edge  of  the  mixing  layer  by  considering  the  entrainment  of 
fluid  into  the  mixing  layer.  Very  near  the  edge  of  the 
turbulent  mixing  layer  the  production  and  dissipation  of 
turbulent  kinetic  energy  is  small  compared  to  convection  and 
diffusion.  Therefore,  the  turbulent  kinetic  energy  equation 
reduces  to 


uifji  +  vifFi  +  p?l?(rkva!i'2)  -  o  <”) 

By  assuming  local  similarity  such  that 


where  a  is  a  constant,  Equation  75  can  be  expressed  as 

h*"  ■  va  5'2>  06) 

Near  the  outer  edge  of  the  mixing  zone  the  change  of 
(Un/a  -  V) ,  which  is  the  entrainment,  with  respect  to  n  can 
be  assumed  to  be  small  compared  to  the  change  of  q' 2  with 
respect  to  n  [24],  Thus,  Equation  76  is  rewritten  as 

nkl(£  o  -  v)5,2l  -  !^<nk  va  5,2>  (76a) 

Integrating  this  equation  gives 
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v ,  =  5-  U  -  V  (77) 

a  o 

The  mass  entrainment  rate  into  one  edge  of  the  mixing 
layer  is 


dQ  = 
dx 


(2irr)k  pUdr 


(78) 


Also /  the  mass  entrainment  rate  into  the  mixing  layer  is 


dQ  .. 
dx 


(27rb)k  p  VD 


(79) 


where  VD  is  the  lateral  velocity  of  the  mass  being  entrained 
into  the  mixing  zone.  Combining  Equations  78  and  79  and 
with  the  aid  of  Leibnitz  rule 


„  T7  db 
D  “  Ub  dx 


U 

W 


Pi,  b'  - 
Kb  o 


8  pUr 
3x 


dr 


and  with  the  aid  of  the  continuity  equation  this  reduces  to 


v  =  u  —  -  V 
D  b  dx  b 


(80) 


where  and  are  the  boundary  values  of  the  velocity  com¬ 
ponents.  Since  similarity  has  been  assumed 

db  _  n 
dx  c 


Then  Equations  77  and  80  give 
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V 


d 


V 


D 


Thus,  the  rate  of  turbulent  diffusion  at  the  boundary  of  the 
mixing  zone  is  equal  to  the  rate  of  mass  entrainment.  This 
observation  permits  use  of  mass  entrainment  measurements  to 
establish  the  turbulent  diffusion  velocity  at  the  boundary. 

In  the  first  regime  of  symmetric  jet  mixing,  mass  is 
entrained  through  the  inner  and  outer  boundaries .  In  the 
second  regime  of  a  symmetric  jet  mass  is  entrained  only 
through  the  outer  boundary  and  the  net  flux  of'  turbulent 
kinetic  energy  being  diffused  across  the  line  of  symmetry  is 
zero . 

In  the  symmetric  jet  the  turbulent  diffusion  velocity 
is  finite  on  the  inner  boundary  in  the  first  regime,  but  as 

previously  noted  it  is  zero  on  the  centerline  in  the  second 

\ 

regime.  Hence,  the  turbulent  diffusion  velocity  at  the  end 
of  the  first  regime  experiences  a  finite  discontinuity,  or  a 
transition  zone  exists  over  which  the  diffusion  is  altered. 
In  either  case  it  is  obvious  that  the  same  function  for  the 
diffusion  velocity  could  not  apply  to  first-regime  mixing  as 
well  as  to  second-regime  mixing. 

Following  Bradshaw  et  al.  [24]  the  turbulent  diffu¬ 
sion  velocity  is  normalized  based  on  the  local  maximum 
absolute  shear  stress  and  a  characteristic  velocity,  U^, 
such  that 


U 


^max 

2 

ch 


(81) 
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The  subscripts,  1  and  2,  denote  the  first  and  second  regimes, 

respectively.  For  the  equilibrium  boundary-layer  problem 

Bradshaw  found  the  diffusion  parameter,  f,  to  be  a  universal 

function.  Choice  of  the  characteristic  velocity  is  subject 

to  some  conjecture  for  the  symmetric  jet,  but  it  is  chosen 

as  (U.  -  U  ),  the  local  maximum  velocity  difference  across 
1  o  J 

the  j  et . 

In  establishing  the  function,  f  ( (r  -  r^)/b) ,  the 
boundary  values  are  established  from  empirical  entrainment 
data.  For  the  high  velocity  mixing  boundary,  data  from  the 
coaxial  jet  experiments  by  Chriss  [18]  and  Paulk  [19]  are 
plotted  in  Figure  3  and  compared  with  Bradshaw's  [24] 
empirical  result  for  the  equilibrium  boundary  layer. 

The  agreement  is  reasonable  and  this  suggests  that  a 
similar  function  exists  for  the  high  velocity  mixing 
boundary ,  i . e . , 


V. 


U 


ch 


=  10 


'max 


hv 


U 


ch 


On  the  inner  boundary  in  the  first  regime  f^(0)  is  -10.  The 
sign  is~negative  since  is  negative  on  the  inner  boundary. 

A  similar  correlation  is  also  suggested  in  Figure  3 
for  the  low  velocity  boundary,  i.e., 


For  U. 
i 


U 


ch 


=  2.5 


T 


max 


LV 


U 


ch 


>  U  ,  f ■  (1)  is  2.5. 
o  1 


With  the  boundary  conditions 
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Figure  3. 


Entrainment  velocity  correlation. 
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established  for  f the  complete  function  was  established  by 
trial  and  error- by  comparing  calculated  and  experimental 
velocity  and  shear  stress  profiles.  The  resulting  diffusion 
function  is  given  in  Figure  4.  Ultimately,  the  diffusion 
function  should  be  established  experimentally.  Negligible 
axial  pressure  gradients  existed  in  the  experiments  and  no 
pressure  gradient  was  assumed  in  the  trial  and  error  cal¬ 
culations  . 

In  the  second  regime  f 2  ( 0 )  is  zero  by  symmetry.  At 
the  outer  boundary'f 2 (1)  was  chosen  as  2.5  to  be  continuous 
with  its  value  in  the  first  regime.  Additional  information 
can  be  derived  from  the  turbulent  kinetic  energy  equation 
for  f2(0).  Consider  the  turbulent  kinetic  energy  equation 
written  in  the  following  form: 


where 


T  =  2.  isLL  +  v  Sail  +  i_ 

2  Sx  2  Sr  rk  3r 


V 


M" 


V. 

4  s,! 


(F1  -  F2>If  +  E<<3'2>3/2  +  D 


(83) 


On  the  centerline  SU/Sr  and  t  are  zero.  Therefore,  T  is  not 
only  zero  on  the  centerline,  but 


ST 

Sr 


=  0 


(84) 
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Figure  4.  Diffusion  function  for  axisymmetric  jet  mixing. 
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which  is  obtained  using  L’ Hospital's  rule.  Realizing  that 
q' 2  and  U  are  symmetric  functions,  Equation  84  reduces 
further  to 


3  V,  q 


i  2 


3r' 


3r  ‘ 


0 


Using  Equation  81  this  gives 


3  2f 
3r  2 


£ 


o 


(85) 


Using  the  boundary  values  established  for  the 
functional  form  was  established  by  trial  and  error.  The 
resulting  function  is  given  in  Figure  4. 

The  negative  values  of  the  function  for  the  first 
regime  represent  inward  diffusion  of  turbulent  kinetic 
energy,'  and  the  positive  values  represent  outward  diffusion 
of  turbulent  kinetic  energy.  In  the  first  regime  the  diffu¬ 
sion  function  is  zero  at  approximately  the  peak  shear  stress 
point. 


III.  DISSIPATION  LENGTH  PARAMETER 

The  dissipation  length  parameter,  L,  in  the  dissipa¬ 
tion  term  of  the  turbulent  kinetic  energy  equation  is 
similar  to  the  dissipation  length  used  by  Bradshaw  et  al. 
[24]  and  discussed  by  Townsend  [26] .  The  parameter  L  is 
also  analogous  to  the  Prandtl  mixing  length  as  can  be  seen 
by  neglecting  the  convection,  diffusion,  and  compressibility 
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terms  in  the  turbulent  kinetic  energy  equation;  hence,  from 
Equation  31 


3/2 


which  is  rearranged  to 


This  shows  that  the  Prandtl  mixing  length  theory  actually 
equates  the  rate  of  local  production  and  local  dissipation 
of  turbulent  kinetic  energy  and  consequently  does  not  incor¬ 
porate  any  history  of  the  turbulent  process.  Including  the 
convection  and  diffusion  terms  in  the  turbulent  kinetic 
energy  equation  does  add  history  to  the  turbulent  structure 
behavior . 

Following  the  approach  of  Lee  and  Harsha  [30] ,  L  is 
taken  to  be  proportional  to  the  width,  b,  of  the  mixing 
layer  such  that 


L  =  |  (86) 

i 

where  K  is  a  constant.  From  Sami’s  data  [20  and  21]  for  an 
incompressible  jet  the  value  of  K  is  plotted  in  Figure  5 
versus  the  nondimen sional  jet  width.  Although  there  appears 
to  be  a  slight  variation  in  K  with  r,  a  constant  and  some* 
what  higher  value  equal  to  0.85  was  found  to  be  more  accept¬ 
able  for  the  calculations  presented  in  Chapter  VI.  In 


46 


1.0 


AEDC-TR-70-134 


Chapter  VI  it  is  shown  that  physically  unrealistic  velocity 
profiles  are  produced  with  K  equal  to  0.6  such  that  the  peak 
velocity  profile  did  not  occur  on  the  centerline.  The 
sensitivity  of  the  mixing  calculation  to  K  will  be  discussed 
in  Chapter  VI. 

It  is  also  of  interest  to  note  how  L  compares  with 

the  dissipation  length  reported  by  Bradshaw  et  al.  [24]  for 

the  boundary  layer.  This  is  shown  in  Figure  6.  Bradshaw's 

dissipation  length  parameter,  denoted  as  L1 ,  is  equivalent 
3/2 

to  L  |a^|  '  .  As  Figure  6  shows,  the  two  functions  are 

remarkably  similar  except  at  the  outer  boundaries.  Bradshaw 
included  the  effects  of  intermittency  which  is  not  included 
in  the  treatment  of  the  mixing  problem. 
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CHAPTER  V 

NUMERICAL  PROCEDURES 

I .  GENERAL  APPROACH 

The  method  of  characteristics  was  used  in  solving  the 
basic  flow  equations.  The  characteristic  slopes  and  corre¬ 
sponding  compatibility  equations  were  numerically  solved  in 
finite  difference  form  on  an  IBM  360-50  computer.  The 
numerical  equations  were  programmed  in  Fortran  IV.  The 
numerical  method  used  for  the  method  of  characteristics  is 
basically  Method  II  of  Ralston  and  Wilf  [39]  with  an  implicit 
numerical  integration  scheme.  This  method  marches  all  the 
characteristics  a  constant  Ax  before  the  next  Ax  increment 
is  taken.  This  makes  the  method  especially  convenient  for 
comparison  of  experimental  data  with  the  calculations  since 
profiles  of  experimental  data  are  generally  taken  and  pre¬ 
sented  at  a  given  x  station. 

II.  SUBROUTINES 


Starting  Line  Subroutine 

The  vertical  velocity  component,  V,  does  not  appear 
in  the  differential  equation  along  the  inclined  character¬ 
istics  (Equation  66) ,  but  only  in  the  differential  equation 
along  the  vertical  characteristic  (Equation  65).  Therefore, 
V  is  calculated  by  direct  integration  along  the  vertical 
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characteristic  with  U  and  t  given  on  the  initial  profile. 


Field  Point  Subroutine 

The  field  point  subroutine  is  schematically  illus¬ 
trated  in  Figure  7.  The  objective  is  to  calculate  U,  V,  and 
t  on  the  new  profile  given  their  values  on  the  old  profile. 
In  symmetric  flow  the  most  logical  scheme  is  to  proceed  with 
the  calculations  from  the  inner  boundary  or  line  of  symmetry 
to  the  outer  boundary.  In  plane  two-dimensional  unsymmetric 
free  shear  flow  the  inner  and  outer  boundary  conditions  are 
related  to  one  another  as  given  by  Equation  73.  In  this 
case  the  calculations  can  be  started  from  either  boundary. 
The  numerical  procedure  is  to  guess  either  or  Vq  and 
iterate*  until  Equation  73  is  satisfied.. 

The  numerical  form  of  the  characteristic  and  com¬ 
patibility  equations  is  expressed  in  Equations  87  through 
91.  The  left  characteristic  equation  is 


r3  -  r2  " 


Ax 
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23 
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V23  + 
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+  K^1  ‘  F2U)23’ 


1/2 


(87) 
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Old  Profile  New  Profile 


Figure  7.  Schematic  representation  of  the  mesh  point 
numerical  method. 
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The  right  characteristic  equation  is 


r4  “  r3  = 
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(88) 


The  compatibility  equations  along  the  characteristics  in 
numerical  form  are  as  follows: 

The  left  compatibility  equation  is 
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where 
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The  right  compatibility  equation  is 
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where 
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and 
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The  vertical  characteristic  equation  is 
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The  double  subscript  notation  indicates  average  values.  For 
example,  U23  is  (U2  +  U3)/2.  The  single  subscript  notation 
is  the  value  of  the  parameters  at  the  designated  mesh  point. 
The  calculation  method  proceeds  as  follows: 

1.  The  position,  r3,  is  first  approximated  by  pro¬ 
jecting  a  left  characteristic  through  point  two  to  point 
three.  The  slope  of  this  characteristic  is  calculated  based 
on  the  average  of  values  at  point  two  and  guessed  values  at 
point  three.  The  first  estimate  of  all  the  parameters  at 
point  three  is  taken  as  the  average  at  the  previously  cal¬ 
culated  point  above  point  two  and  at  point  two. 

2.  A  right  characteristic  is  projected  back  through 
point  three  to  make  the  first  approximation  of  the  position 
of  point  four.  For  the  approximation  of  r^  the  same  average 
values  of  the  parameters  are  used  as  were  used  for  the 
initial*- guess  on  the  characteristic  from  point  two  to  three. 

.3.  By  linear  interpolation,  U^,  V^,  and  are 
determined . 

4.  Using  the  right  and  left  compatibility  equations 
(Equations  89  and  90) ,  the  first  calculated  values  of  U3  and 
r3  are  obtained. 
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5.  Using  the  vertical  compatibility  equation 

(Equation  91),  is  calculated.  Then  the  subroutine 

iterates  through  steps  two  to  five  until  the  calculations 

converge.  The  convergence  criteria  was  arbitrarily  set  at 
-2 

10  per  cent  on  V^,  x-j,  and  r^. 

6.  The  process  is  repeated  at  other  mesh  points  to 
obtain  U,  V,  and  x  over  the  whole  range  of  r. 

7.  The  new  profile  is  then  used  as  input  for  the 
next  profile. 

Inner  Boundary  Subroutine 

For  the  first-regime  case  the  inner  boundary  mesh 

points  are  shown  schematically  in  Figure  8a.  The  right 

characteristic  on  the  boundary  is  line  BD  and  is  arbitrarily 

set  at  99  per  cent  of  the  velocity  difference  between  the 

inner  and  outer  streams  for  U .  >  U  .  The  location  of  the 

3  ° 

boundary  would  be  at  1  per  cent  of  the  velocity  difference 

if  U.  < U  .  With  this  criteron  only  r.  and  x.  need  to  be 
jo  11 

calculated  on  the  new  profile  since  and  U2  are  equal  and 
is  taken  as  zero  in  the  symmetric  flow  problem.  As  was 
mentioned  in  Chapter  III,  is  related  to  Vq  in  plane  two- 
dimensional  unsymmetric  flow  and  an  iteration  is  required  to 
satisfy  this  relation.  The  right  characteristic  equation. 
Equation  88,  and  its  compatibility  equation.  Equation  90, 
are  to  be  simultaneously  solved  for  the  two  unknowns  r^  and 
t^.  At  point  three  the  flow  properties  are  solved  by  the 
field  point  subroutine. 
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Old  Profile  New  Profile 


Old  Profile  New  Profile 


b .  Second-Regime 

Figure  8.  Schematic  of  the  inner  boundary  mesh  points. 
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The  second -regime  inner  boundary  in  symmetric  flow 
requires  an  altogether  different  routine.  The  mesh  points 
are  schematically  shown  in  Figure  8b.  As  was  discussed  in 
Chapter  III,  the  problem  is  to  calculate  only  the  centerline 
velocity  since  the  shear  stress  and  vertical  component  of 
velocity  are  zero  on  the  centerline.  Equation  72  is  used  to 
calculate  the  centerline  velocity,  U  ,  from  conditions 
known  on  the  old  profile.  This  equation  in  finite  differ¬ 
ence  form  is 


U  (U  -  U  ) 
c12  C1  c2 


-( 


i  - 1 


=12 


=12 [dPl 

— w12 


+  (1  +  k)  Ax  (92) 

r  23 

Again  at  point  three  the  flow  properties  are  calculated  by 

the  field  point  routine.  Calculation  of  the  centerline 

velocity,  U  ,  is  an  implicit  calculation  to  be  consistent 
C1 

with  the  regular  field  point  routine.  This  requires  average 
values  of  t  and  r  in  Equation  92,  i.e.,  t ^  an<*  r23' 
respectively.  Therefore,  the  centerline  velocity  calcula¬ 
tion  was  made  in  an  iterative  loop  with  the  calculation  of 
the  properties  at  point  three  until  it  converged. 

Outer  Boundary  Subroutine 

The  outer  boundary  subroutine  is  the  same  in  both  the 
first  and  second  regimes  for  symmetric  flow  and  also  can  be 
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used  for  the  plane  two-dimensional  unsynunetric  flow  case. 

The  mesh  points  of  the  outer  boundary  are  schematically 
shown  in  Figure  9.  In  this  case  the  outer  boundary  is 
arbitrarily  set  at  1  per  cent  of  the  maximum  velocity  differ- 
ence  for  >  UQ.  The  objective  of  this  subroutine  is  to 
calculate  U,  V,  t,  and  r  at  point  three.  The  left  character¬ 
istic  equation  (Equation  87),  its  compatibility  equation 
(Equation  89),  the  vertical  compatibility  equation  (Equation 
91)  and  the  defined  boundary  condition  for  U  are  used  to 
solve  for  the  four  unknowns.  Again  an  implicit  calculation 
procedure  is  used.  The  first  guessed  values  at  point  three 
are  taken  as  those  of  point  two. 

The  procedure  for  stopping  the  field  point  subroutine 
calculations  and  calling  the  outer  boundary  subroutine  is  as 
follows:  If  point  one  on  the  new  profile  is  located  by  the 

field  point  subroutine  such  that  the  right  characteristic 
passing  through  this  point  intersects  the  old  profile  at 
r^  >  point  one  is  discarded  and  point  one  prime  is  used 

in  the  outer  boundary  calculations . 

III.  LOGIC 

The  computer  program  logic  is  illustrated  in  Figure 
10  for  the  symmetric  two-dimensional  jet  mixing  program. 

IV.  STABILITY 

The  Courant-Friedrichs-Levy  criterion 
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jS 


Figure  10.  Flow  diagram. 
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for  a  linear  hyperbolic  system  of  equations  has  so  far 
proved  sufficient  for  numerical  stability.  The  term  Ar  is 
the  step  size  in  the  r-difrection,  Ax  is  the  step  size  in  the 
x-direction  and  (dr/dx)  is  the  maximum  slope  of  either 
the  right  or  left  characteristic.  The  step  size  in  the  com¬ 
puter  program  was  chosen  such  that  this  criterion  was 

\ 

satisfied.  However ,  test  calculations  were  made  for  ft  equal 
to  1.5  and  no  stability  problems  were  encountered.  Bradshaw, 
et  al.  [24]  reports  occasional  instability  near  the  wall  for 
the  boundary-layer  problem  unless  ft  <  0.9.  This  occurred 
when  they  used  a  least-squares  linear  extrapolation  of  the 
shear  stress  to  the  wall. 

Early  in  the  checkout  phase  of  this  computer  program, 
a  fourth-order  Lagrangian  curve  fit  was  used  for  interpo¬ 
lation  on  the  old  profiles.  This  created  an  instability  in 
the  profiles  near  the  centerline.  The  instability  could  be 
recognized  by  the  growth  of  bumps  in  the  profiles,  particu¬ 
larly  the  velocity  profile.  The  problem  was  eliminated  by 
using  a  second-order  Lagrangian  (linear)  curve  fit. 

V .  CONVERGENCE 

A  convergence  criterion  of  0.01  per  cent  was  used  on 
U,  V,  t,  and  r  in  the  iteration  which  proved  sufficient  to 
the  fourth  significant  figure.  With  this  convergence  cri¬ 
terion  an  average  of  four  iterations  for  each  mesh  point  on 
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the  new  profile  was  generally  required. 

VI .  COMPUTATIONAL  TIME 

The  computational  time  depended  on  the  number  of  mesh 
points  chosen  to  start  the  calculation,  the  convergence 
criteria,  and  Ax.  The  number  of  mesh  points  chosen  to  start 
the  calculations  on  the  initial  profile  varied  from  20  to 
40.  In  first-regime  calculations  an  additional  point  is 
picked  up  at  each  new  profile  along  the  inner  boundary.  In 
the  second  regime  the  number  of  mesh  points  remain  constant. 
Generally,  it  was  found  better  to  use  about  twice  as  many 
mesh  points  in  the  inner  half  of  the  mixing  zone  as  in  the 
outer  half  to  obtain  an  adequate  plot  of  the  shear  stress 

r 

profile.  The  program  calculates  along  left  characteristics 
which  monotonically  marches  the  mesh  points  away  from  the 
centerline.  The  program  may  require  adding  more  mesh  points 
near  the  centerline  should  the  distance  between  the  center- 
line  and  first  mesh  point  become  too  great.'  The  calculations 
reported  in  Chapter  VI  did  not  exceed  18  jet  diameters  from 
the  initial  plane  of  mixing,  and  did  not  require  added  mesh 
points . 

The  Ax  used  in  the  mixing  calculations  was  0.12  jet 
diameters  which  satisfied  the  Courant-Friedrichs-Levy 
stability  criterion  in  all  cases  calculated.  This  constant 
value  of  Ax  was  chosen  for  the  convenience  of  comparing  the 
analytical  calculation  with  the  experimental  data  at 
particular  x-stations.  Should  a  variable  Ax  be  chosen  by 
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letting  always  equal  unity  the  program  would  march  faster. 

As  an  example  of  the  computing  time,  about  1.5 
minutes  are  required  for  the  program  to  march  one  jet 
diameter  with  40  mesh  points.  No  attempts  were  made  to 
decrease  the  computing  time  since  no  more  than  15  minutes 
were  generally  required  to  calculate  the  flow  fields  of 
interest. 


VII.  ACCURACY 

The  numerical  accuracy  of  the  computer  program  was 
checked  on  each  calculated  profile  by  inserting  the  results 
into  the  integral  form  of  the  mass  continuity  and  momentum 
equation.  The  details  are  given  in  Appendix  C.  Mass  con¬ 
tinuity  was  always  satisfied  within  2  per  cent  or  better  and 
momentum  was  always  satisfied  within  4  per  cent  or  better. 
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CHAPTER  VI 

THEORETICAL  CALCULATIONS  FOR  THE 
AX I SYMMETRIC  JET 

I.  SELECTION' OF  INITIAL  PROFILE  DATA 

Analytical  calculations  are  made  and  compared  with 
experimental  data  for  hydrogen-air  and  air-air  coaxial  con¬ 
stant  pressure  flow  systems  using  the  method  of  character¬ 
istics  discussed  in  Chapters  III  and  V.  The  hydrogen-air 
experimental  data  is  taken  from  Chriss1  experiments  [18] 
and  the  air-air  experimental  data  is  taken  from  Paulk's 
experiments  [19].  Both  sets  of  data  were  based  on  measure¬ 
ments  made  in  the  same  experimental  apparatus .  Mixing  cal¬ 
culations  are  presented  using  initial  velocity  and  turbulent 
shear  stress  profiles  in  both  the  first  and  second  regimes. 
However,  no  calculations  are  made  and  presented  for  marching 
the  calculations  from  the  first  regime  into  the  second 
regime  because  the  transition  problem  discussed  in  Chapter 
IV  remains  to  be  solved.  The  initial  shear  stress  data 
profiles  were  provided  by  Chriss  and  Paulk  from  their  mean 
flow  measurements  [18  and  19]  using  the  momentum  integral 
method  described  by  Peters  et  al.  [42]. 

For  first-regime  calculations  the  initial  starting 
plane  was  chosen  approximately  2.5  nozzle  diameters  from  the 
initial  mixing  plane,  a  location  where  reasonably  accurate 
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velocity  and  turbulent  shear  stress  profiles  were  obtained. 

The  analytical  calculations  to  be  presented  for 
second-regime  mixing  generally  used  the  second  profile  of 
experimental  data  in  the  second  regime.  The  experimental 
turbulent  shear  stress  profile  was  feared  to  be  inaccurate 
on  the  first  profile  of  experimental  data.  The  turbulent 
shear  stress  data  is  based  on  evaluating  X-derivatives  of  the 
momentum  flux  integral.  In  or  near  the  transition  regime 
these  X-derivatives  are  probably  not  very  accurate  because 
of  the  rapid  change  in  character  of  the  flow  from  the  first 
regime  to  the  second  regime. 

The  analytical  calculations  are  presented  first  for 
hydrogen-air  coaxial  mixing  and  then  for  air-air  coaxial 
mixing. 


II.  HYDROGEN-AIR  COAXIAL  MIXING 

Theoretical  calculations  are  made  and  compared  with 
experimental  data  for  the  case  of  a  central  jet  of  hydrogen 
coaxially  exhausting  with  air  into  a  plenum  at  constant 
pressure.  The  calculations  are  made  for  hydrogen-to-air 
velocity  ratios  ranging  from  2.4  to  6.3.  The  maximum  jet 
Mach  number  is  about  0.6.  Table  I  lists  the  bulk  flow 
properties  for  the  cases  calculated  and  compared  with  experi¬ 
ment.  Variations  in  the  parameters  a^,  K,  and  f  are  made  to 
observe  their  sensitivity  on  the  mixing  calculations  since 
basic  measurements  of  these  parameters  are  not  available 
from  the  experiments.  These  calculations  are  made  for  a 
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TABLE  I 


HYDROGEN-AIR  AXISYMMETRIC 
MIXING  FLOW  CONDITIONS 


V“o 

U  . 
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H./H 
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BTU/lb 

P  -/p 
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1900 

11.8 

1.82  X  103 

0.08 

4.4 

3200 

11.6 

1.78 

0.09 

cn 
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Ul 

3300 

11.9 

1.84 

0.09 

4.6 

2880 

7.7 

1.98 

0.14 
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second-regime  case  since  their  effect  can  be  more  readily 
observed  over  greater  axial  length.  First,  however,  an 
example  of  a  first-regime  hydrogen-air  mixing  calculation  is 
presented. 


First  Regime 

This  example  is  presented  for  a  central  jet-to-outer 
stream  velocity  ratio  of  4.6.  The  calculations  start  at  an 
x/D  of  2.7.  Theoretical  and  experimental  velocity  and 
'turbulent  shear  stress  profiles  are  shown  in  Figures  11  and 
12  at  3; 57  and  4.53  nozzle  diameters  from  the  jet  efflux. 

The  diffusion  function,  f^,  used  in  these  calculations 
is  that  shown  in  Figure  4,  page  44.  The  parameter  K  is 
chosen  equal  to  0.85.  Two  analytical  functions  for  the 
parameter  a^  are  used  in  these  calculations  for  comparison. 

In  one  case  it  is  assumed  equal  to  -0.17  and  in  the  other 
case  it  is  a  variable  given  by 


0.9 


+  0.1, 


Z  <  Z 
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=  1.0 


1  ’m 


z  > 
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(94) 


where  |a^lm  and  Zq  are  equal  to  0.17  and  0.4,  respectively. 

This  function  for  a1  is  a  modified  form  of  Equation 
74.  Equation  74  could  not  be  used  for  first-regime  mixing 
since  the  inner  boundary  is  defined  where  the  shear  stress 
is  finite;  consequently,  the  turbulent  kinetic  energy  is 
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a.  x/D  =3.57 

Figure  11.  First-regime  velocity  profiles  for  hydrogen-air 
mixing . 
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b.  x/D  =4.53 
Figure  11.  (continued) 
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b.  x/D  =4.53 
Figure  12.  (continued) 
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finite.  Therefore,  Equation  94  is  used  as  an  estimated 
approximation  for  the  parameter  a^.  The  two  functions  used 
for  a^  do  not  produce  enough  difference  in  the  velocity 
profiles  to  be  seen  on  a  plot  of  this  scale.  The  shear 
stress  profile  near  the  inner  boundary  of  the  mixing  zone  is 
more  sensitive  to  a^,  particularly  near  the  end  of  the  first 
regime.  The  proper  shape  of  the  shear  stress  profile  is  not 
produced  by  the  analytical  calculations  near  the  inner 
boundary  for  either  function  of  a1  although  the  function 
given  by  Equation  94  seems  to  give  the  better  result.  The 
problem  is  not  unexpected  since  the  transition  region  is 
being  approached.  Also,  a  better  function  for  a^  is  probably 
needed  since  Equation  94  is  only  an  estimated  approximation. 
Some  improvement  in  the  diffusion  function  may  also  be 
needed.  Producing  the  proper  shear  profile  at  the, end  of 
the  first  regime  is  necessary  before  the  marching  integra¬ 
tion  can  proceed  from  the  first  into  the  second  regime.  An 
improper  shear  stress  profile  at  the  end  of  the  first  regime 
will  result  in  an  improper  decay  of  the  centerline  velocity 
in  the  second  regime.  With  the  exception  of  the  problem 
just  discussed,  the  velocity  and  shear  stress  profiles  cal¬ 
culated  agree  reasonably  well  with  the  experimental  data. 

Second  Regime 

Second-regime  calculations  are  presented  for  the  flow 
conditions  given  in  Table  I,  page  68.  First,  the  sensitivity 
of  the  calculations  of  velocity  and  shear  stress  profiles 
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and  the  centerline  velocity  and  the  peak  shear,  stress  decay 
to  the  parameters  a^,  K,  and  fj  is  presented  in  Figures  13 
through  29.  The  experimental  data  for  the  hydrogen-air 
velocity  ratio  of  4.4  is  used  for  comparison  and  the  initial 
profiles  were  selected  at  6.5  jet  diameters  from  the  initial 
plane  of  mixing.  For  the  best  values  of  a^,  K,  and  f 2  found 
in  these  calculations,  calculations  for  the  other  velocity 
ratio  cases  will  be  presented. 

The  empirical  expression  for  a^  used  in  the  calcula¬ 
tions  is  the  analytical  function  given  by  Equation  74  in 
Chapter  IV.  The  peak  level  of  a^  is  established  by  | a.^ | m 
and  lateral  position  of  this  peak  is  established  by  ZQ.  The 
sensitivity  of  a1  on  the  calculations  will  be  observed  by 
choosing  different  constant  values  of  la. I  and  Z  while  K 
and  f2  are  held  constant.  In  Figures  13  and  14  the  center- 
line  velocity  and  peak  shear  stress  decay  is  observed  to  be 
relatively  insensitive  to  The  choice  of  |a1|m  equal 

to  0.17  based  on  Sami's  data  (discussed  in  Chapter  IV)  for 
an  incompressible  air  jet  is  quite  satisfactory  in  this 
case.  This  implies  that  the  relation  between  the  turbulent 
shear  stress  and  turbulent  kinetic  energy  is  essentially 
independent  of  density  effects,  because  there  is  an  order- 
of-magnitude  variation  in  density  across  the  two  coaxial 
flows.  From  this  result  perhaps  one  also  could  expect 
little  or  no  compressibility  effects  on  the  function  a^. 
Bradshaw  and  Ferriss  [25]  found  this  to  be  the  case  for  the 
adiabatic  compressible  boundary  layer. 
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AEDC-TR-70-1 34 


Figure  15.  Effect  of  la]_lm  on  the  velocity  profile 
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Figure  16.  Effect  of  la]_lm  °n  the  shear  stress  profile. 
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Figure  18.  Effect  of  ZQ  on  the  peak  shear  stress. 
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Figure  23.  Effect  of  K  on  the  velocity  profile 
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O  Experiment  (Ref.  18) 
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Figure  26.  Effect  of  f,  on  the 


J _ I _ 1 _ I _ I 

10  12  14  16  18 


> 

m 

o 


shear  stress 


AEDC-TR-70-134 


90 


AEDC-TR-70-134 


Figure  28.  Effect  of  f2  on  the  shear  stress  profile. 
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The  effect  of  la^lm  on  the  velocity  and  shear  stress 

profiles  at  an  x/D  of  14.5  is  shown  in  Figures  15  and  16, 

respectively.  The  velocity  profile  is  more  sensitive  to 

la. I  near  the  centerline,  whereas  the  shear  stress  is  more 
'I'm 

sensitive  away  from  the  centerline.  At  this  axial  position 
the  peak  shear  stress  is  overpredicted  about  30  per  cent. 

No  conclusive  reason  can  be  given  for  this  although  there 
are  several  factors  which  could  be  the  cause.  For  instance, 
the  theory  incorporates  the  assumption  of  unity  Prandtl  and 
Schmidt  number  which  Chriss  [18]  has  observed  to  be  less 
than  unity  for  hydrogen-air  mixing  .j  A  slight  axial  variation 
of  |a.l  ,  ZQ/  and  K  may  be  closer  to  physical  reality. 

The  sensitivity  of  the  centerline  velocity  and  peak 
shear  stress  calculations  to  ZQ  is  shown  in  Figures  17  and 
18.  A  relatively  weak  effect  is  observed  here  but  as  seen 
in  Figure  19  the  velocity  profile  shape  is  strongly  affected 
near  the  centerline.  In  fact,  a  value  of  0.3  for  ZQ  is 
physically  incorrect  for  this  case  since  the  peak  velocity 
should  occur  on  the  centerline.  Obviously,  an  improper 
balance  and  rate  of  change  of  convection,  diffusion,  pro¬ 
duction,  and  dissipation  of  turbulent  kinetic  energy  near 
the  centerline  is  created  in  the  calculations.  A  reasonable 
velocity  profile  shape  is  produced  with  ZQ  equal  to  0.4  or 
0.5.  The  better  choice  is'  0.4.  This  gives  reasonable 
agreement  with  experimental  data  and  also  fits  Sami's  data 
very  well  in  the  use  of  Equation  74.  This  value  of  ZQ  is 
not  always  the  best  choice  for  other  mixing  flows  as  will  be 
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shown  later  in  this  chapter.  In  Figure  19,  page  82,  varia¬ 
tion  in  Zq  is  observed  to  produce  a  slight  change  in  the 
peak  shear  stress. 

Analytical  calculations  are  made  for  K  equal  to  0.6, 
0.85,  and  1.0.  The  effect  of  this  variation  in  K  on  the 
centerline  velocity  and  peak  shear  stress  decay  is  shown  in 
Figures  21  and  22,  pages  84  and  85.  The  effect  on  the 
lateral  velocity  and  shear  stress  profiles  is  shown  in 
Figures  23  and  24,  pages  86  and  87.  Too  much  velocity  decay 
is  produced  near  the  centerline  for  K  equal  to  0.6  as 
observed  in  Figure  23.  As  with  a^,  the  velocity  profile  is 
most  sensitive  to  K  near  the  centerline.  Of  the  three  cases 
calculated  K  equal  to  0.85  gives  the  best  agreement  with 
experimental  data.  This  is  a  somewhat  higher  value  of  K 
than  that  calculated  from  Sami's  experimental  data  and  pre¬ 
sented  “in  Figure  5,  page  47. 

f 

The  diffusion  velocity  function,  f2,  was  formulated 
following  the  approach  of  Bradshaw  et  al.  [24}  as  discussed 
in  Chapter  IV.  The  sensitivity  of  the  velocity  and  turbu¬ 
lent  shear  stress  calculations  is  obtained  by  scaling  f^  by 
three  different  constants  such  that  fjd)  equals  1.5, 

2.5,  or  3.5.  The  effect  on  the  centerline  velocity  and  peak 
shear  stress  decay  is  shown  in  Figures  25  and  26,  pages  88 
and  89.  The  velocity  profile  is  mainly  affected  near  the 
centerline  as  observed  in  Figure  27,  page  90.  This  center- 
line  effect  is  caused  by  the  influence  that  the  diffusion 
velocity  has  on  the  mixing-layer  thickness  which  in  turn 
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affects  and  the  energy  dissipation  term.  Both  of  the 
latter  parameters  are  functions  of  the  mixing-layer  thick¬ 
ness.  In  Figure  28,  page  91,  the  shear  stress  profile  is 
observed  to  be  only  slightly  sensitive  to  a  variation  of  f 2* 
A  good  measure  of  the  proper  choice  of  the  diffusion  func¬ 
tion  near  the  mixing  boundary  is  its  prediction  of  the 
mixing-layer  growth  rate.  The  results  of  these  calculations 
are  shown  in  Figure  29,  page  92.  Two  conclusions  are  drawn. 
First,  the  form  of  the  diffusion  function  chosen  (Equation 
81)  predicts  the  mixing-layer  growth  rate  very  well. 

i 

Secondly,  the  function,  f 2,  proposed  in  Figure  4,  page  44, 
produces  reasonable  velocity  and  shear  stress  profiles.  For 
f2(l)  equal  to  3.5  the  mixing  boundary  spreads  too  fast  and 
for  f 2  (1)  equal  to  1.5  it  spreads  too  slowly. 

To  summarize  the  sensitivity  calculations  analytical 

and  experimental  velocity  profiles  are  presented  in  Figure 

30  at  8.6,  12.5,  and  14.5  jet  diameters  from  the  initial 

plane  of  mixing.  The  choices  of  la. I  ,  Z  ,  K,  and  f_(l) 

i  m  o  z 

are  0.17,  0.4,  0.85,  and  2.5,  respectively.  Although  the 
calculated  velocity  profiles  agree  well  with  the  experi¬ 
mental  data  with  this  combination  of  values  of  the  param¬ 
eters,  other  better  combinations  of  these  parameters 
may  exist.  A  detailed  experimental  study  of  the  turbulent 
structure  is  needed  to  better  assess  these  parameters. 

Calculations  on  the  remaining  velocity  ratio  cases 
listed  in  Table  I,  page  68,  are  presented  in  Figures  31 
through  33.  The  "best"  values  of  K,  and  found 
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Figure  30.  Lateral  velocity  profiles. 
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Figure  31.  (continued) 
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c.  Velocity  Profile 
Figure  31.  (continued) 
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d.  Turbulent  Shear  Stress  Profile 
Figure  31.  (continued) 
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b.  Peak  Turbulent  Shear  Stress  Decay 
Figure  32.  (continued) 
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b.  Peak  Turbulent  Shear  Stress  Decay 
Figure  33.  (continued) 
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in  the  sensitivity  calculations  are  used  and  are  found  to 
work  very  well.  However,  the  best  value  for  ZQ  is  observed 
to  vary  from  one  case  to  the  other.  This  suggests  that  a^ 
is  not  a  universal  function  of  Z.  An  estimate  for  ZQ  equal 
to  Z  at  t  was  found  to  work  occasionally  but  not  always. 
Still,  a  better  analytical  function  for  a^  is  desired  near 
the  line  of  symmetry.  A  choice  of  a^  equal  to  -0.17  over 
the  outer  mixing  zone  works  well. 

III.  AIR-AIR  COAXIAL  MIXING 

The  air-air  mixing  calculations  are  made  for  a  central 
jet-to-outer  stream  velocity  ratio  of  8.0.  The  nominal 

“\ 

inviscid  jet  and  outer  stream  velocities  are  400  and  50 
ft/sec,  respectively.  The  maximum  Mach  number  in  the  mixing 
flow  field  is  approximately  0.35,  low  enough  to  be  con¬ 
sidered  incompressible  flow.  Both  first-regime  and  second- 
regime  calculations  are  presented  and  compared  with  experi¬ 
mental  data.  The  functions  for  the  parameters  la-jJ^,  K, 
and  f2  used  for  the  hydrogen-air  mixing  calculations  are 
used  here. 

First  Regime 

The  initial  profiles  of  experimental  data  were  taken 
at  2.6  jet  nozzle  diameters  from  the  nozzle  exit.  Velocity 
and  shear  stress  profiles  are  calculated  and  compared  with 
experimental  data  at  4.4  jet  nozzle  diameters  from  the 
initial  plane  of  .mixing,  a  position  less  than  one-half 
nozzle  diameter  from  the  end  of  the  first  regime.  These 
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results  are  shown  in  Figures  34  and  35.  As  with  the 
hydrogen-air  first-regime  calculations  the  same  two  func¬ 
tions  for  a^  are  used  in  these  analytical  calculations. 
Similar  to  the  hydrogen-air  results  no  detectable  differ¬ 
ences  are  observed  in  the  velocity  profiles  for  the  two 
functions  of  a.^  used.  The  shear  stress  profile  also  exhibits 
an  improper  shape  near  the  inner  boundary.  Otherwise,  the 
theoretical  prof ile  calculations  agree  well  with  the  experi¬ 
mental  velocity  and  shear  stress  profiles.  These  observa¬ 
tions  generally  agree  with  the  results  of  the  hydrogen-air 
mixing  calculations. 

Second  Regime 

One  second-regime  theoretical  calculation  is  made  for 
coaxial  air  streams.  The  experimental  data  for  a  jet-to- 
outer  stream  velocity  ratio  of  8.0  is  taken  from  Reference 
19.  Plots  of  the  centerline  velocity  and  peak  shear  stress 
decay  are  shown  in  Figures  36  and  37.  The  theoretical  cal¬ 
culations  were  made  for  la, I  and  K  equal  to  0.17  and  0.85, 

1  i  m 

respectively,  the  same  numerical  values  as  were  used  for  the 
hydrogen-air  calculations.  The  diffusion  function,  is 

also  the  same  as  that  used  for  the  hydrogen-air  calculations. 
The  centerline  velocity  agrees  well  with  the  experimental 
data.  The  peak  shear  stress  experimental  data  are  scattered; 
consequently,  good  agreement  of  the  theory  with  experimental 
data  is  not  expected.  The  results  do  indicate  that  the  peak 
shear  stress  decay  is  properly  predicted. 
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Figure  35.  First-regime  turbulent  shear  stress  profile  for 
air-air  mixing. 
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IV.  OTHER  THEORETICAL  CALCULATIONS 


Compressibility 

In  the  development  of  the  turbulent  kinetic  energy 
equation  term  D  (Equation  28)  appears  as  the  result  of 
including  density  fluctuation  terms.  If  this  term  and  the 
other  density  fluctuation  terms  are  neglected  the  form  of 
the  turbulent  kinetic  energy  equation  left  is  the  incom¬ 
pressible  form.  Calculations  were  made  for  hydrogen-air  and 
air-air  mixing  to  determine  the  effect  of  D  on  the  results. 
The  effect  is  very  small.  Differences  only  in  the  fourth 
significant  figure  in  the  calculations  of  U,  V,  and  t  were 
observed.  Therefore,  the  term  can  be  neglected  at  least  for 
subsonic  incompressible  hydrogen-air  and  air-air  mixing. 

For  supersonic  mixing  it  may  be  necessary  to  include  D  in 
the  calculations. 

Vertical  Component  of  Velocity 

A  typical  result  of  the  calculated  vertical  component 
of  velocity  is  shown  in  Figure  38  for  hydrogen-air  and  air- 
air  coaxial  mixing.  In  the  air-air  mixing  case  the  stream¬ 
lines  have  positive  slope  near  the  centerline  and  negative 
slope  in  the  outer  mixing  zone.  Therefore,  the  mass  flux  is 
decreasing  near  the  line  of  symmetry  although  the  net  mass 
flux  is  increasing  in  the  mixing  zone.  In  the  hydrogen-air 
mixing  case  a  very  different  behavior  is  observed.  The 
streamlines  have  negative  slope  over  most  of  the  mixing 
zone.  The  mass  flux  near  the  centerline  is  increasing  even 
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Figure  38.  Calculated  vertical  component  of  velocity 
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though  the  velocity  decays.  This  is  primarily  caused  by  the 
displacement  of  the  light  hydrogen  with  the  heavier  air 
through  the  action  of  turbulent  mixing. 
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CHAPTER  VII 

THEORETICAL  CALCULATIONS  FOR  THE 
TWO-DIMENSIONAL  SYMMETRIC  WAKE 

The  results  of  a  brief  study  on  extending  the  theo¬ 
retical  method,  presented  in  Chapters  II  and  III,  to  a  two- 
dimensional  symmetric  wake  is  presented  in  this  chapter. 
Theoretical  calculations  are  made  and  compared  with  the 
experimental  data  of  Reference  44  for  an  incompressible 
turbulent  wake  behind  a  thin  flat  plate.  Before  presenting 
the  results,  a  brief  discussion  is  given  on  the  choice  of 
functions  for  the  turbulent  kinetic  energy-shear  stress 
parameter,  the  diffusion  velocity,  and  the  energy  dissipa¬ 
tion. 


I.  PARAMETER  MODELING 

Basically,  the  same  functional  forms  as  were  developed 
for  the  coaxial  jet  mixing  problem  are  used  for  the  wake 
problem.  I 

Turbulent  Kinetic  Energy -Shear  Stress  Parameter 

The  turbulent  kinetic  energy -shear  stress  parameter, 
a1#  is  in  part  arbitrarily  chosen  to  be  given  by  Equation  74 
with  |a1!m  and  ZQ  equal  to  0.15  and  0.25,  respectively.  The 
value  of  ZQ  selected  is  strictly  arbitrary.  The  value  of 
la. I  selected  is  based  on  that  recommended  by  Reference  40 
for  plane  two-dimensional  wakes. 
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The  diffusion  velocity  parameter  used  is  given  by 

Equation  81  where  the  characteristic  velocity,  U^,  is 

chosen  as  U  -  U  .  The  diffusion  function,  f,  used  in  the 
o  c 

calculations  is  arbitrarily  selected  to  be  the  same  as  given 
in  Figure  4,  page  44,  for  f ^  but  scaled  up  by  a  factor  of 
four  to  provide  the  proper  high-velocity  mixing  boundary 
condition,  i.e.,  to  give  f ^  (1)  a  value  of  10.0. 

Dissipation  Length  Parameter 

The  dissipation  length  parameter,  L,  is  again  given 
by  Equation  86.  From  trial-and -error  calculations  a 
constant  value  of  0.53  for  K  was  found  to  be  satisfactory 
for  the  wake. 


II.  RESULTS  AND  DISCUSSION 

The  initial  starting  profile  for  the  wake  calcula¬ 
tions  is  taken  at  the  end  of  the  flat  plate.  The  general 
numerical  procedures  are  the  same  as  are  described  in 
Chapter  V  for  second-regime  mixing  calculations  with  one 
exception.  Once  the  position,  y,  of  the  left  characteristic 
nearest  the  line  of  symmetry  exceeds  its  initial  value  by  a 
factor  of  two,  a  new  mesh  point  was  added  by  linear  inter¬ 
polation  halfway  between  the  line  of  symmetry  and  the  first 
mesh  point.  The  results  of  these  calculations  are  presented 
in  Figures  39  through  42.  The  centerline  velocity  is 
plotted  versus  the  distance  from  the  end  of  the  flat  plate. 
This  distance  is  nondimensionalized  by  the  turbulent  boundary 
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Figure  40.  Maximum  turbulent  shear  stress  in  the  wake  of  a  thin  flat  plate. 
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Figure  41.  Velocity  profile  in  the  wake  of 
a  thin  flat  plate. 
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Figure  42.  Turbulent  shear  stress  profile  in  the  wake  of  a 
thin  flat  plate. 
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layer  momentum  thickness,  0q,  at  the  end  of  the  flat  plate 
calculated  by  Reference  [44]  to  be  0.53  cm.  The  theory  pre¬ 
dicts  the  centerline  velocity  increase  very  well.  The 
maximum  turbulent  shear  stress  decay  is  not  predicted  very 
well  near  the  end  of  the  flat  plate  but  improves  farther 
downstream.  The  velocity  profile  at  414  momentum  thick¬ 
nesses  downstream  of  the  flat  plate  is  in  reasonable  agree¬ 
ment  with  the  experimental  data.  The  calculated  shear  stress 
profile  is  also  in  reasonable  agreement  with  experimental 
data  with  the  maximum  deviation  occurring  in  the  outer 
mixing  2one. 

These  results  show  that  the  theoretical  method  is 
applicable  to  wake  calculations.  However,  the  best  func¬ 
tions  for  a^  K,  and  f  remain  to  be  determined.  Refinement 
of- these  parameters  based  on  experiments  is  necessary  before 
completely  satisfactory  results  can  be  expected. 

The  computer  program  was  also  modified  to  calculate 
the  mixing  flow  field  for  an  unsymmetric  two-dimensional 
plane  jet  or  plane  wake.  This  program  includes  an  iteration 
on- the  vertical  component  of  velocity  to  satisfy  the  condi¬ 
tion  required  by  Equation  73 .  So  far ,  numerical  diffi¬ 
culties  have  been  encountered  in  trying  to  calculate  flows 
where  the  turbulent  shear  changes  sign,  such  as  occurs  in  an 
unsymmetric  wake  behind  an  airfoil. 
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CHAPTER  VIII 

CONCLUDING  REMARKS  AND  RECOMMENDATIONS 

I.  CONCLUDING  REMARKS 

The  results  of  the  calculations  presented  in  Chapters 
VI  and  VII  show  that  the  turbulent  kinetic  energy  equation 
transformed  into  a  transport  equation  for  turbulent  shear 
stress  is  potentially  useful  for  calculating  inhomogeneous 
turbulent  two-stream  mixing  flows  and  wake  flows.  By  using 
the  turbulent  kinetic  energy  equation  the  history  of  the 
turbulent  structure  is  included  in  the  calculations.  The 
eddy  viscosity  and  Prandtl  mixing  length  theories  are  in¬ 
adequate  in  this  respect  because  history  is  excluded.  The 
key  to  successfully  using  the  turbulent  kinetic  energy  equa¬ 
tion  lies  in  properly  modeling  the  diffusion  and  dissipation 
terms  and  obtaining  a  realistic  relation  between  the  turbu¬ 
lent  shear  stress  and  turbulent  kinetic  energy.  The  con¬ 
vective  or  flux  model  used  in  these  studies  for  the  diffusion 
of  turbulent  kinetic  energy  works  well.  A  different  diffu¬ 
sion  function  is  required  in  the  first  and  second  regimes  of 

i 

axisymmetric  mixing  because  of  the  basic  difference  in  the 
character  of  the  two  flow  fields.  The  same  diffusion  func¬ 
tions  work  equally  well  for  hydrogen-air  and  air-air  mixing, 
thereby  suggesting  that  the  functions  are  universal  in  their 
respective  regimes  of  application. 
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The  model  for  the  dissipation  function  used  in  this 
study  is  analogous  to  the  dissipation  expression  derived  for 
isotropic  turbulence.  By  setting  the  dissipation  length,  L, 
proportional  to  the  mixing  layer  width,  the  model  chosen 
for  the  dissipation  function  works  very  well.  The  parameter, 
K,  equal  to  0.85  served  to  calculate  hydrogen-air  and  air- 
air  coaxial  mixing  flows  equally  well.  These  results  imply 
that  the  parameter  L  used  in  this  study  exhibits  a  universal 
behavior.  A  value  of  K  equal  to  0.53  worked  better  for  the 
two-dimensional  symmetric  wake. 

The  turbulent  shear  stress  was  assumed  to  be  linearly 
proportional  to  the  turbulent  kinetic  energy.  The  parameter 
a.^  is  a  variable  function  for  axisymmetric  jet  mixing.  On 
the  line  of  symmetry  its  value  is  zero.  A  constant  value  of 
0.17  for  | a^^  over  the  outer  region  of  the  mixing  zone 
works  well  for  jet  mixing.  In  the  inner  mixing  region  |a^| 
ranges  from  zero  to  0.17.  A  better  formulation  of  a^  in 
this  regrion  is  desired  opposed  to  that  used  in  this  study. 

In  formulating  the  theory  the  assumptions  of  unity 
Prandtl  and  Schmidt  numbers  were  made.  The  velocity  and 
shear  stress  profiles  and  centerline  velocity  and  peak  shear 
stress  decay  agree  well  with  experimental  data  for  jet  mixing 
even  though  experimental  results  show  that  the  Prandtl  and 
Schmid fnumbers  are  not  unity. 

The  turbulent  kinetic  energy  approach  presented  in 
this  study  is  not  yet  a  useful  engineering  tool  for  calcu¬ 
lating  mixing  flow  fields.  Since  the  method  is  designed  to 
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include  flow  history  in  the  calculations  initial  profile 
data  are  needed  to  begin  the  calculations.  Experimental 
data  have  been  used  in  this  study,  but  analytically  derived 
starting  profiles  are  needed  for  engineering  computations.  , 

A  study  on  the  feasibility  of  using  analytically  derived 
starting'  profiles  is  desired. 

Potentially  the  turbulent  kinetic  energy  approach  is 
a  method  for  studying  the  effects  of  initial  disturbances  in 
turbulent  flow  on  mixing  phenomena.  One  can  envision  use  of. 
the  method  to  study  the  effect  of  initial  boundary  layers  on 
turbulent  wakes  and  to  study  methods  for  retarding  or 
enhancing  the  mixing  in  turbulent  flows. 

II.  RECOMMENDATIONS  FOR  FUTURE  STUDY 

Additional  experimental  confirmation  of  the  turbulent 
kinetic  energy  approach  used  in  this  study  is  desired. 
Studies  are  recommended  for  both  jets  and  wakes  to  compare 
the  convection,  diffusion,  production,  and  dissipation  of 
turbulent  kinetic  energy  measured  experimentally  with  that 
calculated  by  the  analytical  models  used  in  this  study. 

Such  a  study  would  be  useful  in  further  substantiating  these 
analytical  models  and  perhaps  would  suggest  improved  analyt¬ 
ical  models.  In  particular,  a  more  satisfactory  analytical 
model  for  the  parameter  a1  is  needed  for  jet  mixing. 

Studies  should  continue  toward  solving  the  transition 
regime  between  the  first  and  second  regimes  in  symmetric 
two-stream  mixing.  The  proper  analytical  models  for  a^  and 
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the  diffusion  are  needed  before  success  can  be  expected. 

Including  nonunity  Prandtl  and  Schmidt  numbers  in  the 
theoretical  method  is  desired  to  permit  proper  calculation 
of  the  energy  and  species  concentration  profiles.  This  can 
be  accomplished  by  simultaneously  solving  the  parabolic  form 
of  the  energy  and  species  continuity  equations  along  with 
the  global  continuity,  momentum,  and  turbulent  kinetic  energy 
equations.  A  finite  difference  grid  procedure  would  be 
required  opposed  to  the  method  of  characteristics  because 
this  system  of  equations  is  not  hyperbolic.  Another  rather 
involved  approach  would  be  to  use  the  transport  equations , 
analogous  to  the  turbulent  kinetic  energy  equation,  for  the 
fluctuating  energy  and  the  species.  These  equations,  how¬ 
ever,  contain  six  additional  parameters  [43]  which  need  to 
be  experimentally  determined  and  appropriately  modeled. 

Application  of  the  theory  to  supersonic  mixing  should 
be  studied.  Compressibility  is  already  included  in  the 
basic  flow  equations.  The  main  question  is  whether  or  not 
the  parameter  modeling  can  be  equally  applicable  to  incom¬ 
pressible  and  compressible  flow. 

Finally  it  is  recommended  that  mixing  flows  with 
pressure  gradient  be  studied  using  the  turbulent  kinetic 
energy  approach. 
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APPENDIX  A 


ORDER  ANALYSIS  OF  THE  CONTINUITY ,  MOMENTUM,  AND 
TURBULENT  KINETIC  ENERGY  EQUATIONS 

To  facilitate  simplification  of  the  equations  of 
motion  and  the  turbulent  kinetic  energy  equation  relative 
estimates  of  the  derivatives  and  of  the  perturbation  com¬ 
ponents  are  needed.  If  L  is  chosen  as  a  length  scale  con¬ 
forming  to  the  variation  of  parameters  in  the  X-direction 
and  similarly  H  for  the  Y-direction  and  it  is  assumed  that 
£/L  < 1  (Order  of  0.1),  then 


3 


9y 


Also  Ug  is  chosen  as  a  characteristic  velocity  such  that  it 
represents  the  velocity  difference  across  the  mixing  zone, 
and  U/Ug  is  assumed  to  be  of  order  unity.  It  is  further 
assumed  that 


3r 

3y 
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Then 
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and 


u*  -  O'  (u„/I7l) 

s 


(A5) 


I.  CONTINUITY 


From  Equations  17  and  AS 

I7l)  (A6) 

The  order  of  V  is  found  from  the  continuity  equation.  First 
note  that 


PU 


(A7) 


then  retaining  the  highest  order  terms,  the  continuity  equa¬ 
tion  is  expressed  as 


|£E  +  !_(pV  +  yW)  -  0 


(A8) 


From  Equations  A5  and  A6 
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fpO.l 


which  is  of  the  same  order  as  3pU/3x.  This  gives 


V 


(A9) 


Note  that  p 'V'  and  pV  are  of  the  same  order  which  requires 
that  both  terms  must  be  retained  in  the  continuity  equation. 
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II.  MOMENTUM 

Assuming  that  the  viscous  stresses  in  the  momentum 
equation  (Equation  7,  Chapter  II)  are  small  compared  to  the 
turbulent  stresses,  the  X-  and  Y-components  of  momentum  are 
written  with  the  order  of  each  term  written  below.  The  bars 
over  the  mean  flow  parameters  have  been  dropped  for  con¬ 
venience  . 
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Y-coraponent 
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If  all  terms  of  an  order  less  than  pU2/L  are  neglected  the 

S 

X-  and  Y-components  of  the  momentum  equation  reduce  to 
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Y-component 


3P  =  3pV' 2 
3y  3y 

The  latter  equation  can  be  integrated  directly  to  give 

P  =  P  -  pV1"2 
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(All) 


where 
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This  term  is  neglected  since  it  is  an  order  less  than  the 
order  of  the  terms  in  the  X-component  of  momentum. 


III.  TURBULENT  KINETIC  ENERGY 


Before  an  order  analysis  is  performed  on  the  turbulent 
kinetic  energy  equation  (Equation  10,  Chapter  II)  the  con¬ 
vection  term  (I)  and  normal  stress  term  (III)  are  combined 
and  rearranged  with  the  aid  of  the  continuity  equation  as 
follows : 
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With  this  rearrangement  the  turbulent  kinetic  energy  equa¬ 
tion  is  restated  as 


(f)  (g) 


Each  of  the  first  five  terms  are  expanded  in  two-dimensional 
Cartesian  coordinates  with  the  order  of  the  respective  terms 
listed . 
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kLJ  #  L 

lJ 

3/2 


1 


+  (pUTVT  +  p'V’U')^  +  (pvr2  +  yryr2)|X 


pUgfoi  2 


fel 2  !^fAl5/2  pUsfM  PUsfiLl3/2 
lLJ  '  L  W  ~[Lj  '  IT [l 


3  .  3 

IT'ILJ '  ~L~\L\ 
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Term  (e) 


P'u'  +  J  pq'2U'  +|  p*q' ZU' j 


p0s(r) 2 '  pDs(e]3/2'  pus(r) 


3 

3y 


Fv1 


1 

1 


pq 1  2V  ’  + 


k 


T\P 


1 

Jt 


2 

9 


pU 


s 


(4) 


2 

; 


Term  (f) 


+  P 


,  3  V ' 
3y 


pUsfA] 

L  [LJ 


Term  (g ) 

In  reviewing  the  order  of  each  of  the  preceding 

terms,  the  highest  order  term  occurs  only  in  term  d. 

Therefore,  the  viscous  dissipation  must  be  of  the 

same  order  and  should  be  retained. 

Retaining  all  terms  of  the  order  of  puH/L2  or  higher 

s 

reduces  the  turbulent  kinetic  energy  equation  to 


pU  (pV  +  p'V1  )  qTZ/2 

2  3x  2  3y  3y 

+  §|  +  V  ^  +  PuTZ  +  (piPv1'  +  p'U'V'jjp- 
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+  pv'  2  +  UFvr  +  i  pq rzVT  +  4  p',q,iV 


9y  9y 


9  a 


+  pT|YI  +  u!  pa=  o 

3y  l  9Xj 


(A12) 


If  only  terms  of  the  order  of  (puVL)  (A/ L)  or  higher  are 

S  ' 

retained  the  turbulent  kinetic  equation  reduces  to 


(pinr  +  p'u'v1)^  + 


i  ^  9  a !". 

5  pq^J  +  “55^  =  «  <A13> 


These  three  terms  remaining  are  the  production,  diffusion, 

i 

and  dissipation  of  turbulent  kinetic  energy.  The  convection 
term  is  at  least  half  of  an  order  of  magnitude  smaller  than 
the  diffusion  which  itself  is  half  an  order  of  magnitude 
smaller  than  production  and  dissipation.  The  pTrUTVT  term  is 
included  in  the  production  term  because  it  is  of  the  same 
order  as  the  diffusion.  This  causes  no  problem  with  the 
solution  to  the  equations  because  (pU1 V'  +  p 'U'V' )  can  be 
treated  as  one  variable. 

Since  the  ensuing  objective  of  this  study  will  be  to 
include  as  much  physics  as  practical  into  the  solution  of 
the  turbulent  jet,  the  form  of  the  turbulent  kinetic  energy 
equation  given  in  Equation  A12  will  be  used  as  the  starting 
point.  The  normal  stresses  in  term  (d)  and  term  (f)  will  be 
neglected.  Terms  (b)  and  (c)  will  be  defined  as 


TTrT 


pD  =  p'U 


v 

v  ayj 


l 

7 


9prVr  q1 


9y 


(A14) 
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The  form  of  the  turbulent  kinetic  energy  equation  remaining 
is 


k  pu  +  |(pv  +  +  (purvr  +  P'U'V')|^ 


9y 


_  1  _  i  _ 

P '  V '  +  J  pq'  ZV 1  +  i  p 'q ' ZV ' 

•  i 


- 5F7T 

+  Ui  -3^  +  "D  ’  0 


(A15) 


This  form  of  the  turbulent  kinetic  energy  equation  is  used 
in  this  study  of  the  inhomogeneous  two-stream  jet  mixing 


problem. 


p.. 


~  O 


142 


AEDC-TR-70-134 


APPENDIX  B 

DEVELOPMENT  OF  THE  CHARACTERISTIC  AND 
COMPATIBILITY  EQUATIONS 


Using  vector  notation  and  subscript  notation  for  the 
differentials,  consider  the  quasilinear  hyperbolic  differ¬ 
ential  equation 


AW  +  BW  +  C  =  0 
x  r 


where  A  and  B  are  square  matrix  coefficients  and  W  and  C  are 
vectors.  By  definition  of  the  differential 

dW  =  W  dx  +  W  dr 
x  r 

Define  A  to  be  a  characteristic  direction  of  the  system  of 
equations,  (Equation  Bl) ,  such  that 


,  _  dx 
A  "  dF 


(B3) 


patibility  equations  are  found  by  combining  Equations  Bl  and 
B2  to  give 


B  C  -  <B\  -  A)  Wx 


(B4) 


Let  D  be  a  vector  (yet  to  be  defined) ,  and  form  the  inner 
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product  of  both  sides  of  Equation  B4  with  D, 

(B  If  +  C)  *  D  =  <BX  -  A)  wx  *  D 

=  W  ■  (BX  -  A)T  D  (B5) 

X 

T 

where  (BX  -  A)  is  the  transpose  of  (BX  -  A) .  If  D  can  be 
chosen  such  that  the  right-hand  side  of  Equation  B5  is 
identically  zero  with  respect  to  W  ,  then  Equation  B5  will 
contain  derivatives  only  in  the  characteristic  direction,  X, 
and  therefore,  will  be  the  compatibility  equation  in  that 
direction.  Such  a  D  does  indeed  exist,  and  is  obtained  as  a 
non-trivial  solution  of  the  linear  algebraic  system, 

(BX  -  A)T  D  =  0  (B6 ) 


Non-trivial  solutions  of  Equation  B6  are  guaranteed  as  a 
result  of  Equation  B3.  The  technique  outlined  is  used  to 
derive  the  characteristic  directions  and  compatibility 
equations  from  the  continuity,  momentum,  and  turbulent 
kinetic  energy  equations.  Denote  the  unknown  vector,  W,  as 


W  = 


r  \ 

u 

V 

T 


(B7) 


From  the  partial  differential  equations  (Equations  54,  55, 
and  56) ,  the  matrices  A  and  B  and  vector  C  are 


A  = 


fl  +  U(F1  -  F2)  0 

U  0 

0  0 


0 

u 


2a 


1' 
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B  = 


V(F1  -  F2) 

V  -  T  (F1  -  F2) 


^<Pl  ‘  F2>  -  1 


-  T 


1 

0 


0 

-1 


V  +  V, 


tU 


2a, 


2a, h 


fl  -  F„U ; 

l  J-  J-  J 

T 

2 

2a  i 

h  J 

(B8 ) 


X  =  0 


(B9) 


or 


dr 

dx 


(B9a) 


Corresponding  to  this  solution  the  unknown  vector  D  is 
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1  +  U(F,  - 


f2)  0 


Dj^tl  +  -  F2)  ]  +  D2D  - 


217  D3  ’  0 


A  non-trivial  solution  is 


D1  =  -U/[l  +  U(F1  -  F2)] 


°2  =  1 


d3  =  0 


(BIO) 


Substituting  this  solution  into 


dr 


+  C  •  D  =  0 


(Bll) 


gives  the  vertical  compatibility  equation.  To  illustrate 


V<F1  -  P2) 

V  -  T  (F1  -  F,2) 


T  2^(F1  -  F2J  -  1 


[  T  ]  1  -  F2u 

2a!  E 


V  +  V 


d  tU  dx 


2a,  2a,  h  dr 

■i-  J- 
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Then 


Thus 


b^  +  c  = 

dr 


[V(F  -  F  ^  +  ^+  Zk  +  Ii^l 
V(F1  F2}dr  +  dr  +  r  +  P  dxj 


[V  -  T  (F.  -  F0)  y. 

1  1  2  dr  dr  r 


(y  -  V 

h  dP' 

V  \ 

P  dx. 

2i7(Fi  -  V 


-  X 


1  f  T  ^ 

_  n-  \  T  1 

•l  -  F2U] 

dU 

i — 1 
(0 
CM 

h  J 

dr 

-  T 


1  "  F2U 

l2aJ 

h  J 

§!•  +  G 
dr 


E  ^  +  C 
dr 


U 

[1  +  U(FX  "  F2)] 
+  [V  “  T  (F^  -  F^)] 


* 

V(F1  - 

. 

dU  _  dr 
dr  ""  dr 


F  )  —  + 
2;dr 


-  -  k  + 
r 


dV 

dr 


U  dP 
P  dx 


h  dP 
P  dx 


=  0  (B12) 

which  reduces  to  Equation  61. 

The  compatihility  relations.  Equations  62,  correspond¬ 
ing  to  nonzero  X  are  obtained  in  a  similar  manner. 
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APPENDIX  C 

CONSERVATION  OF  MASS  AND  MOMENTUM 

Conservation  of  mass  and  momentum  are  checked  to 
determine  the  accuracy  of  the  method  of  characteristics  used 
in  the  calculation  of  the  mixing  flow  field.  The  check 
method  uses  the  integral  form  of  the  continuity  and  momentum 
equations.  First,  consider  the  mass  continuity  check. 

I.  MASS  CONTINUITY 


The  mass  continuity  equation  is  integrated  between 
the  inner  and  outer  boundaries 


(Cl) 


The  second  term  is  integrated  directly  to  give 


9pUrk 

3x 


dr  + 


„  k 
p  V  r 
o  o  o 


p.V.r?  =  0 
l  1 


(C2) 


Using  Leibnitz's  rule  the  first  term  is  rearranged  as 


r 

/ 


3pUr  Hr- 

r  dx  r 

i  i 


r 

=  a-/° 
dx  J 


pUrkdr 


-  I  PUr 


v  dr 
k  o 


o  o  o  dx 


.  dr. 

p.U.rk  -s— - 
l  l  dx 


(C3 ) 
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Combining  Equations  C2  and  C3 


a-f‘ 

dx  J 

r . 

1 


pUrkdr 


.  dr  ,  dr. 

rT  k  o  „  k  1 
p  U  r  j —  -  p.U.r.  - — 
po  o  o  dx  Ki  1  1  dx 


-  (p  V  rk  -  p. V.rk) 
o  c  i  i ' 


(C4 ) 


Equation  C4  is  integrated  with  respect  to  x,  giving 


■ 

r . 
i 


pUrkdr 


L 

-r 

x  r. 


pUrkdr 


X 

■/ 

x_ 


k  dro  k  dri 
p  U  r  -  p.U.r.  ^ 
Fo  o  o  dx  Ki  i  x  dx 


dx 


X 

-I 


(p  V  rk 
o  o  o 


p . V . r .  )  dx 

ill 


x . 

j 


(C  5 ) 


where  xg  is  the  position  of  the  initial  profile.  Letting  P 
be  defined  as 


•L 

■/ 


R  S  ]  pUrkdr 

q  1 


r . 
1 


JL 

/' 

ri 


pUrkdr 


x 

/ 


n  k  dro  _  „  n  k  dri 
poUoro  dx  piUiri  dx 


dx 


X  X 

o  s 
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X 


-/ 

x_ 


"oVo 


piViri*  dx 


(C6) 


Rg  is  evaluated  at  each  axial  station  by  substituting  into 

the  integrals  the  values  of  the  parameters  calculated  by  the 

method  of  characteristics.  By  letting  x  always  be  the 

position  of  the  initial  profile,  the  variation  of  R  from 

PI 

unity  will  indicate  the  error  accumulated  by  the  method  of 
characteristics . 


II .  MOMENTUM 

The  differential  forms  of  the  mass  continuity  and 
momentum  equations  (Equations  33  and  34)  are  combined  to 
give 


3pU2rk  3pVUrk  _  3rkT  _  dP 

3x  3r  3r  dx 


(C7) 


By  integrating  across  the  mixing  zone 


3  pU 


2  k 


3x 


dr  +  p  V  U  rk  -  p . V . U . rk 
Ko  o  o  o  1111 


_k  k  (  .  dP 

roTo  riTi  ro  ri  dx 


(C8) 


The  first  term  is  integrated  by  Leibnitz's  rule. 


X 

I 

r . 

i 


o  artTT2  k 


=  —  f' 

dx  J 


pU2rkdr  - 


,  dr 

tt  k  o 
p  U  r  j — 
o  o  o  dx 


r . 
i 
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.  dr. 
*iuiri  dST 


(C9) 


By  combining  Equations  C8  and  C9, 


r 

r 

dx  J 


,a-k 


pUrdr  = 


i  dr  ,  dr . 

p  u2rk  _  p  ,u?r^  —=■ 
Ho  o  o  dx  1  1  dx 


r . 
i 


+  rkT  -  rkx  -  — (r  -  r  ) 
+  roTo  riTi  dxiro  riJ 


-  p  V  U  rk  +  p.V.U.rk 
Ho  O  O  O  1  1  1 


(CIO) 


and  integrating  with  respect  to  x  gives 


r 


pU2rkdr 


r 

l 


,2  k. 


pU^r  dr 


x 


A 

■I 


,  dr 

,,2  k  o 
p  Ur  3 — 
Ko  o  o  dx 


,  dr. 

p.U?rk  3—^1 
i  i  dx  ) 


dx 


/ 


+  1 

O  O 


k.  , 

r . t . ) dx 

l  l 


x 

-/ 


(r  - 


r  . )  dx 
l  dx 


x 


/ 


pUVrdx  + 
o  o  o  o 


x 

/p . U . V . rkdx 
1111 


(Cll) 


By  definition  let 
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XT 

■/' 


R  =  I  pU2rkdr 
m  1 


r 

r 

r . 

1 


pU2rkdr 


dr. 


X=  Xe 

s  s 


/  (poUoro  33T 


.  dr. 

p.U?r.  -r-i- 
l  i  dx 


dx 


X 

/ 


s 


+  1  {roTo 

o  o 


k  .  , 
r . t . ) dx 
l  l 


-/ 


p  U  V  r  dx 
o  o 


s 


X  X 

+  j  p.U.V^dx  -  j  (ro  -  r. )  ||  dx 


(C12) 


Again,  the  integrals  at  each  of  the  appropriate  x  stations 

are  evaluated.  The  accumulated  error  in  momentum  calculated 

by  the  method  of  characteristics  is  determined  by  comparing 

R  to  unity, 
m 

In  the  computer  program  R^  and  R^  are  calculated  at 
each  x  station  so  that  the  loss  or  gain  in  mass  or  momentum 
can  be  observed. 
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